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Abstract 

We  consider  a  phase  field  model  for  the  formulation  and  solution  of  topology  op¬ 
timization  problems  in  the  minimum  compliance  case.  In  this  model,  the  optimal 
topology  is  obtained  as  the  steady  state  of  the  phase  transition  described  by  the 
generalized  Cahn-Hilliard  equation  which  naturally  embeds  the  volume  constraint 
on  the  amount  of  material  available  for  distribution  in  the  design  domain.  We 
reformulate  the  model  as  a  coupled  system  and  we  highlight  the  dependency  of 
the  optimal  topologies  on  dimensionless  parameters;  also,  we  discuss  the  issue 
of  mesh  dependency  of  the  solution.  We  consider  Isogeometric  Analysis  for  the 
spatial  approximation  which  facilitates  encapsulating  the  exactness  of  the  repre¬ 
sentation  of  the  design  domain  in  the  topology  optimization  and  is  particularly 
suitable  for  the  analysis  of  phase  field  problems.  We  demonstrate  the  validity  of 
the  approach  and  numerical  approximation  by  solving  two  and  three-dimensional 
topology  optimization  problems. 

Keywords:  Topology  optimization;  minimum  compliance;  phase  field,  model;  Isogeo¬ 
metric  Analysis. 

1  Introduction 

In  engineering  it  is  often  desired  to  apply  some  optimization  techniques  to  the  design 
of  a  structure,  component  or  device.  Other  than  sizing  [9,  108]  and  shape  optimiza¬ 
tion  techniques  [52,  69,  96],  a  significant  contribution  is  given  by  topology  optimization 
[10,  12,  14,  84,  86],  which  represents  the  fundamental  form  of  optimization;  indeed,  to¬ 
pology  optimization  aims  at  finding  the  optimal  distribution  of  a  material  in  a  design 
domain  such  that  an  objective  functional  is  minimized  under  certain  constraints.  The 
minimum  compliance  case  represents  the  most  common  topology  optimization  problem, 
for  which  the  goal  is  to  generate  the  globally  stiffest  structure  by  distributing  only  a  lim¬ 
ited  amount  of  material  in  the  design  domain  [10,  12];  additionally,  another  interesting 
problem  consists  in  generating  the  lightest  structure  under  stress  constraints,  see  among 
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the  others  e.g.  [23,  36,  75].  Historically,  topology  optimization  has  been  used  principally 
for  structural  static  problems  based  on  a  linear  elastic  model,  but  many  other  cases  have 
also  been  successfully  considered.  For  example,  this  is  the  case  for  applications  in  fluid 
dynamics  [1],  heat  conduction  [45],  vibration  [58],  multiphysics  [92]  and  bioengineering 
[110];  also  topology  optimization  has  been  used  for  shell  structures  [65,  84]  and  with 
different  material  models  as  in  [88]  for  elastoplastic  structures. 

In  most  cases,  topology  optimization  problems  are  defined  in  simplified  geometries, 
typically  rectangles,  representing  the  design  domain.  Even  if  this  is  a  reasonable  starting 
point,  in  many  cases  it  would  be  interesting  to  perform  topology  optimization  for  a 
part  or  component  of  a  structure  for  which  an  initial  design  already  exists.  Since  it  is 
common  practice  in  engineering  to  represent  geometries  with  Computer  Aided  Design 
(CAD)  technologies,  which  are  based  on  NURBS  [80]  or,  more  recently,  T-splines  [89], 
it  is  desirable  to  include  the  exact  representation  of  the  design  domain  in  the  topology 
optimization  procedure.  However,  in  current  practice,  the  numerical  approximation 
scheme  used  for  topology  optimization,  typically  the  Finite  Element  Method  (e.g.  [32, 
54,  83]),  requires  the  approximation  of  the  design  domain  and  disconnects  the  analysis, 
and  hence,  the  optimization  from  its  geometrical  representation. 

In  general,  the  capability  to  embed  the  CAD  geometric  representation  in  the  analy¬ 
sis  and  optimization  provides  not  only  accuracy  advantages,  but  also  has  the  potential 
to  considerably  improve  the  efficiency  of  the  overall  design  procedure.  The  importance 
of  establishing  a  suitable  link  between  the  optimization  and  the  CAD  representation  is 
recognized  and  discussed  in  [15]  for  shape  optimization,  in  particular  for  shells  struc¬ 
tures;  the  authors  propose  a  procedure  combining  design  modeling,  structural  analysis 
and  optimization,  for  which  these  tasks  are  coordinated  by  means  of  a  program  system 
named  Computer  Aided  Research  Analysis  Tool  [16]  made  available  to  the  designer.  In 
[84]  shape  optimization  problems  are  solved  by  considering  the  position  of  the  control 
points  of  B-spline  [80]  as  design  variables  together  with  adaptive  refinement  strategies; 
a  similar  procedure  is  extended  to  topology  optimization  problems,  combining  repeated 
optimization  steps  with  B-spline  approximations  of  the  optimal  topologies  and  adaptive 
refinement.  In  [69]  the  relation  between  CAD  and  shape  parametrization  is  discussed 
for  shape  optimization,  especially  for  fluid  dynamics;  additionally,  in  [109]  manipulation 
of  the  splines  is  used  to  generate  optimal  geometries.  Also,  in  [63]  topology  optimization 
problems  have  been  solved  by  using  control  points  of  B-spline  curves  as  design  variables 
in  an  approach  combing  shape  optimization  and  hole  nucleation. 

Isogeometric  Analysis,  a  generalization  of  Finite  Element  Analysis  for  which  basis 
functions  are  defined  by  NURBS  or  T-splines  [33,  55],  provides  the  possibility  to  embed 
the  exact  CAD  representation  of  the  design  domain  in  topology  optimization,  in  addition 
to  exhibiting  several  other  advantages  [8,  34,  39,  46].  Isogeometric  Analysis  has  already 
been  introduced  successfully  and  discussed  for  shape  optimization  in  [29,  50,  70,  105] 
and  we  believe  that  it  also  represents  a  potentially  effective  numerical  approximation 
method  for  topology  optimization  problems.  Recently,  Isogeometric  Analysis  has  been 
used  in  [91]  to  solve  design  optimization  problems  to  generate  optimal  two-dimensional 
structures  by  means  of  a  procedure  based  on  trimmed  curves;  this  concept  is  further 
extended  in  [90]  for  topology  optimization  problems.  In  this  manner  the  final  optimal 
structure  is  represented  by  NURBS  and  T-splines  and  directly  linked  to  the  CAD  rep¬ 
resentation  without  the  need  of  additional  postprocessing  of  the  topology  optimization 
result.  However,  even  if  this  represents  a  great  advantage  and  ideal  situation,  the  to¬ 
pology  optimization  results  appear  to  be  strongly  dependent  on  the  specific  approach 
used  to  generate  the  trimmed  curves  and  surfaces.  Additionally,  B-spline  bases  are 
considered  in  [61]  for  two-dimensional  topology  optimization  problems.  Ideally,  a  com- 
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prehensive  design  optimization  procedure  based  on  Isogeometric  Analysis  could  be  used 
to  provide  an  optimized  structure  from  an  initial  design  domain  passing  through  topo¬ 
logy  optimization,  geometry  generation  and  shape  optimization,  while  maintaining  the 
centrality  of  the  geometry  in  the  overall  procedure. 

In  its  original  formulation,  topology  optimization  is  a  distributed  and  discrete  valued 
problem  [12],  for  which  only  areas  of  material  and  void  are  allowed  without  intermedi¬ 
ate  states.  However,  this  formulation  leads  to  many  difficulties  both  from  the  analysis 
and  the  numerical  points  of  view,  and  it  requires  efficient  discrete  optimizers,  see  e.g. 
[98].  The  most  popular  approach  to  overcome  this  difficulty  is  based  on  the  material 
distribution  concept,  for  which  the  design  variable  corresponds  to  a  density  function 
smoothly  representing  the  distribution  of  the  material  in  the  design  domain,  with  inter¬ 
mediate  values  between  the  pure  material  and  void  states  allowed.  In  this  framework  a 
possible  approach  is  the  homogenization  method  [4] ,  for  which  the  macroscopic  proper¬ 
ties  of  the  material  are  deduced  from  the  microscopic  properties  of  the  porous  material 
represented  by  the  density  function;  the  first  numerical  approximation  for  an  homog¬ 
enized  material  was  presented  in  [10].  However,  in  this  approach  solutions  appear  to 
have  an  elevated  number  of  microscopic  holes  and  microstructures  which  are  undesired 
from  a  manufacturability  point  of  view,  when  pure  material  and  void  states  are  re¬ 
quired.  In  order  to  obtain  these  kinds  of  optimal  topologies,  the  intermediate  states  can 
be  penalized  by  choosing  suitable  interpolation  schemes  for  the  dependency  of  macro¬ 
scopic  material  properties  on  the  density  function;  in  the  case  of  isotropic  materials, 
the  Solid  Isotropic  Material  Penalization  (SIMP)  model  is  the  most  successfully  used 
[11,  12,  67,  87].  Typically,  topology  optimization  problems  in  this  approach  are  solved 
with  suitable  constrained  optimization  techniques  and  with  low-order  Finite  Element 
approximation  for  the  density  function,  often  piecewise  constant  over  the  elements. 
However,  additional  stabilization  and  filtering  techniques  need  to  be  introduced  at  the 
level  of  the  numerical  approximation  in  order  to  remove  or  reduce  the  well-known  mesh 
dependency  and  checkerboard  phenomena  [12,  59]  which  affect  the  topology  optimiza¬ 
tion  results.  Different  techniques  have  been  considered  to  solve  topology  optimization 
problems,  among  these  are  Evolutionary  Structural  Optimization  (ESO)  [111]  and  Bidi¬ 
rectional  ESO  methods  [113],  heuristic  procedures  based  on  the  identification  of  regions 
of  material  with  high  and  low  contributions  to  the  stiffness  of  the  structure;  also,  other 
optimization  strategies  based  on  the  removal  of  material  by  the  evaluation  of  topological 
derivatives  have  been  adopted  [71].  However,  in  general,  among  the  drawbacks  of  these 
formulations  there  is  the  strong  dependence  of  the  optimal  solution  on  the  particular 
optimizer  utilized  and  its  settings. 

Recently,  the  use  of  the  level  set  method  [74]  has  been  proposed  to  solve  topology 
optimization  problems  [3,  5,  35,  106].  In  this  approach,  the  introduction  of  a  level  set 
function,  which  is  associated  with  the  density  function,  avoids  directly  tracking  the 
boundaries  between  the  material  and  the  void;  the  optimal  solution  is  then  obtained  as 
the  evolution  in  time  of  the  level  set  function,  for  which  an  optimizer  is  no  longer  needed. 
However,  topological  changes  are  uni-directional,  in  the  sense  that  holes  can  only  be 
removed  in  the  design  domain  and  inner  front  creation  requires  additional  numerical 
techniques  [76,  112];  also,  similar  to  other  level  set  methods,  repeated  reinitializations 
of  the  level  set  function  are  required  while  numerically  solving  the  problem. 

An  alternative  approach  to  topology  optimization  is  provided  by  the  multiphase 
formulation,  where  the  distribution  of  two  phases,  representing  the  material  and  void, 
inside  the  design  domain  is  described  by  a  smooth  function  which  coincides  with  the 
density  function.  The  geometrical  information  associated  with  the  optimal  topology  is 
then  deduced  from  the  sharp  interfaces  between  the  two  phases,  which  are  represented  by 
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thin  layers.  The  definition  of  topology  optimization  problems  in  a  multiphase  approach 
has  recently  been  introduced  for  design  dependent  loads  in  [20]  (and  [21]),  for  problems 
with  stress  constraints  in  [23]  and  for  the  minimum  compliance  case  in  [107,  115];  further 
extensions  are  also  considered  in  [102].  The  concept  at  the  basis  of  this  approach  is  that 
the  objective  functional  to  be  minimized  is  penalized  by  means  of  additional  terms 
controlling  the  interfaces  and  the  decomposition  of  the  pure  phases,  which  are  typical  of 
multiphase  problems  [6,  25,  40].  In  particular,  the  introduction  of  the  interface  term  for 
penalization  allows  the  definition  of  a  well-posed  topology  optimization  problem  [23], 
and  it  provides  at  the  discrete  level  optimal  solutions  not  affected  by  mesh  dependency 
and  checkerboard  phenomena.  Still  the  problem  is  formulated  as  an  optimization  one, 
for  which  the  optimal  topology  depends  in  general  on  the  optimizer  used. 

Further,  the  topology  optimization  problem  in  the  multiphase  approach  can  be  trans¬ 
formed  into  a  phase  field,  problem  for  which  the  optimal  topology  is  obtained  as  the 
steady  state  of  the  phase  transition;  at  the  basis  of  this  formulation  there  is  the  reinter¬ 
pretation  of  the  penalized  objective  functional  introduced  for  the  multiphase  approach 
as  a  total  free  energy.  Traditional  phase  field  models  are  represented  by  the  Cahn- 
Hilliard  [25,  26,  27]  and  Cahn-Allen  [6]  equations,  which  have  been  introduced  in  met¬ 
allurgy  to  describe  phase  segregation  in  binary  alloy  systems.  More  recently,  phase  field 
approaches  have  been  successfully  considered  to  provide  mathematical  models  for  prob¬ 
lems  in  different  disciplines;  for  example,  there  are  models  for  crack  propagation  [18,  66], 
also  with  Cahn-Hilliard  equation  [97],  image  segmentation  [104]  and  cancer  and  tumor 
growth  [41,  72].  The  role  of  the  phase  field  approach  for  topology  optimization  consists 
in  obtaining  separated  phases,  material  and  void,  divided  by  thin  and  sharp  interfaces 
for  which  the  distribution  of  the  material  in  the  design  domain  is  determined  by  the 
optimization  considerations.  In  this  sense,  this  resembles  the  case  of  the  Cahn-Hilliard 
equation  with  elastic  misfit,  for  which  the  distribution  of  the  phases  partially  takes  into 
account  the  elastic  properties  of  the  materials;  see  e.g.  [42].  In  topology  optimization, 
this  effect  has  to  assume  a  leading  role,  and  the  distribution  of  the  material  depends 
on  the  objective  function  of  the  topology  optimization  problem.  Phase  field  models  for 
topology  optimization  have  been  considered  firstly  in  [107,  114,  115]  for  the  minimum 
compliance  case,  also  for  multimaterials  problems;  a  nonlinear  fourth-order  generalized 
Cahn-Hilliard  equation  is  derived  and  successfully  solved  for  two  and  three-dimensional 
problems  by  using  a  multigrid  algorithm  [114,  115]  and  with  an  approach  which  par¬ 
tially  decouples  the  phase  and  elasticity  equations.  The  mass  conservation  property 
of  the  Cahn-Hilliard  equation  is  conveniently  used  to  naturally  take  into  account  the 
mass/volume  constraint  associated  with  the  minimum  compliance  problem.  More  re¬ 
cently,  a  similar  approach  based  on  the  Cahn-Allen  equation  has  been  proposed  in  [102] 
for  shape  and  topology  optimization,  even  if  without  the  capability  to  introduce  topo¬ 
logical  changes.  In  [30]  a  model  based  on  the  diffusion-reaction  equation,  with  analogies 
with  the  Cahn-Allen  equation,  is  introduced  for  minimum  compliance  problems  with 
an  augmented  Lagrangian  approach  to  take  into  account  the  mass/volume  constraint. 

Phase  field  models  show  many  similarities  with  the  level  set  approach,  however,  they 
allow  to  naturally  include  hole  nucleation  in  the  formulation  and  avoid  reinitializations 
of  level  set  functions  while  numerically  solving  the  problem. 

In  this  work  we  formulate  the  topology  optimization  problem  for  the  minimum  com¬ 
pliance  case  by  using  a  phase  field  approach  following  [114],  since  we  feel  that  this  kind 
of  formulation  provides  several  advantages.  Namely,  it  has  the  ability  to  naturally  deal 
with  topological  changes,  to  provide  geometrical  information,  and  to  completely  de¬ 
scribe  the  topology  optimization  problem  at  the  continuous  level,  without  the  necessity 
of  introducing  any  ad  hoc  numerical  techniques  at  the  discretization  stage;  also,  since 


4 


the  problem  of  optimization  is  converted  to  a  phase  transition  problem,  the  need  to  use 
an  optimizer,  and  hence  the  dependence  of  the  solution  on  its  settings,  is  eliminated 
and  replaced  with  the  choice  of  a  suitable  time  approximation  scheme.  We  rederive  the 
generalized  Cahn-Hilliard  equation  starting  from  the  SIMP  and  multiphase  approaches 
on  the  basis  of  energy  considerations,  following  the  usual  procedure  for  the  derivation  of 
phase  field  models,  and  we  highlight  the  parametric  dependency  of  the  problem  intro¬ 
duced  by  the  penalization  of  the  objective  functional.  Also,  we  rewrite  the  phase  field 
model  as  a  coupled  system  of  phase  and  elasticity  equations,  we  provide  its  dimension¬ 
less  form,  and  we  characterize  it  in  terms  of  dimensionless  parameters.  We  discuss  the 
choice  of  the  parameters  used  to  penalize  the  objective  functional  and  the  mesh  depen¬ 
dence  of  the  optimal  solution,  and  we  extend  the  continuation  method  [20,  23,  71],  an 
optimization  strategy  based  on  the  sequential  solution  of  optimization  subproblems,  to 
the  phase  field  case. 

For  the  numerical  solution  of  the  topology  optimization  problem  in  the  phase  field 
approach  we  consider  Isogeonretric  Analysis  for  the  spatial  approximation,  since  we  be¬ 
lieve  that  it  provides  several  benefits  for  the  solution  of  this  kind  of  problem  in  analogy 
with  [46,  48]  for  phase  field  problems.  Firstly,  Isogeometric  Analysis  encapsulates  the 
CAD  representation  of  the  design  domain  in  the  topology  optimization  while  providing 
geometric  flexibility;  also,  it  ensures  robustness,  high-order  accuracy,  and  the  capability 
to  easily  use  compactly  supported  high-order  basis  functions  for  the  approximation  of 
the  nonlinear  fourth-order  generalized  Cahn-Hilliard  equation,  whose  numerical  solu¬ 
tion  necessitates  spatially  C 1  -continuous  functions.  For  the  time  approximation  we  use 
the  generalized-a  method  [31]  together  with  a  time-adaptive  scheme  that  allows  a  very 
efficient  solution  of  the  problem,  which,  in  analogy  to  other  phase  field  models,  exhibits 
fast  and  intermittent  variations  in  time.  We  show  the  effectiveness  of  the  proposed  pro¬ 
cedure  by  solving  two  and  three-dimensional  topology  optimization  problems  in  design 
domains  defined  by  NURBS  geometries. 

This  work  is  organized  as  follows.  In  Sec.  2  the  SIMP  and  multiphase  approaches 
for  topology  optimization  in  the  minimum  compliance  case  are  recalled.  In  Sec.  3  the 
derivation  of  the  standard  Cahn-Hilliard  phase  field  model  is  recalled  in  anticipation  of 
the  presentation  in  Sec.  4  of  the  phase  field  model  for  topology  optimization.  In  Sec.  4 
the  generalized  Cahn-Hilliard  equation  is  derived,  the  coupled  system  presented,  and 
a  dimensional  analysis  performed  in  order  to  highlight  the  dependency  of  the  problem 
on  dimensionless  parameters.  In  Sec.  5  we  present  the  numerical  approximation  scheme 
based  on  Isogeometric  Analysis  and  the  generalized-a  method  with  time  adaptivity.  In 
Sec.  6  we  discuss  the  dependency  of  the  optimal  solution  on  the  initial  distribution  of 
material  and  the  choice  of  the  parameters  upon  which  solutions  depend;  also,  we  present 
the  continuation  method  in  the  context  of  the  phase  field  approach.  In  Sec.  7  we  pro¬ 
vide  and  discuss  numerical  results  for  two  and  three-dimensional  topology  optimization 
problems.  Finally,  conclusions  are  presented  in  Sec.  8. 
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Figure  1:  Representation  of  a  two-dimensional  design  domain  Q,  boundaries  Tjj,  T n, 
surface  force  h  and  body  force  f. 


2  Topology  Optimization  in  the  Minimum  Compli¬ 
ance  Case 

In  this  section  we  introduce  the  topology  optimization  problem  in  the  minimum  com¬ 
pliance  case  by  means  of  the  SIMP  and  the  multiphase  approaches,  which  represent  the 
basis  for  the  definition  of  the  phase  transition  model  of  Sec.  4.  Standard  notation  is  used 
through  this  work  to  denote  the  Sobolev  spaces  of  functions  with  Lebesgue  measurable 
derivatives  and  norms;  see  e.g.  [2]. 

2.1  The  SIMP  approach 

Let  us  start  by  introducing  a  material  density  function  p  =  p(x)  to  represent  the  dis¬ 
tribution  of  a  given  material  at  any  generic  point  x  in  a  design  domain  fi  C  of 
dimension  d  =  2,3.  By  convention,  p  =  1  indicates  the  presence  of  the  material,  while 
p  =  0  corresponds  to  regions  of  Q  where  the  material  is  absent,  which  we  will  refer  to 
as  void;  intermediate  states  of  p  between  0  and  1  are  allowed  and  indicate  regions  of 
“soft”  material.  Also,  we  require  that  0  <  p  <  1,  since  values  outside  this  range  do 
not  correspond  to  meaningful  representations  of  the  material  distribution.  We  adopt  a 
formulation  based  on  the  linear  elastic  theory  for  small  displacements  with  an  isotropic 
material  [54],  whose  properties  are  fully  described  by  the  symmetric  elastic  tensor  Co 
which  depends  on  the  Young’s  modulus  Eq,  the  Poisson  ratio  isq  and  the  dimension  d  of 
the  problem  (for  d  =  2,  both  the  plane-stress  and  plane-strain  cases  can  be  considered). 

The  SIMP  approach  is  based  on  the  concept  that  the  properties  of  the  material 
depend  on  the  density  function  p  for  which  the  elastic  tensor  C  (/?)  is  a  function  of  p;  in 
particular,  we  can  write: 

C(p)  =  g(p)C0,  (1) 

with  g{p)  a  suitable  function  introducing  the  homogenization  of  the  elastic  properties 
depending  on  the  distribution  of  the  material  in  fi.  The  function  g{p)  assumes  a  crucial 
role  in  the  definition  of  the  SIMP  method,  since  the  quality  of  the  topology  optimization 
results  strongly  depend  on  it;  we  will  return  on  this  point  later.  It  follows  that  the  stress 
tensor  is  dependent  on  p  as: 

d(p,u)  :=C(p)e(u),  (2) 

where  e(u)  is  the  strain  tensor  associated  to  a  given  displacement  u.  We  observe  that 
the  stress  tensor  d(p,  u)  is  symmetric  and  linearly  dependent  on  u;  the  superscript 
is  used  to  indicate  the  dependency  on  both  p  and  u  as  independent  variables. 
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The  elastic  problem  in  strong  form  consists  in  finding  the  displacement  u,  for  a  given 
material  density  function  p,  such  that: 


-V  •  cr(p,u)  =  f 

in  f 2, 

u  =  0 

on  Td, 

a(p,  u)n  =  h 

on  T at, 

p  given, 

where  To  C  1  is  the  Dirichlet  partition  of  the  design  domain  boundary  dfl  where  the 
displacement  is  imposed,  while  T^r  :=  dCl\Tjj >  is  the  part  of  the  boundary  where  the 
surface  force  h  is  applied  (traction  or  pressure);  for  the  sake  of  simplicity  we  assume  a 
null  displacement  on  Tp.  Also,  f  is  the  body  force  acting  in  the  domain  Q  which  we 
consider  independent  on  p.  An  example  is  presented  in  Fig.  1  for  a  two-dimensional 
design  domain. 

Let  us  introduce,  in  view  of  the  weak  form  of  the  elastic  problem  (3),  the  func¬ 
tion  spaces  V  :=  |v  G  [iL1(0)]d  :  v|rD  =  oj  and  H  :=  {</  €  L°°(fl)  :  g{p)  € 
moreover,  we  define  the  residual  Ru( u; p)(v)  G  R  such  that: 

i?u(u;p)(v)  :=  f  a(p,  u)  :  e(v)  d£l  —  I  f  •  v  dfl  —  ®  h  •  v  cfT  /v ,  (4) 

v  o  </  J  r  jv 

where  we  assume  that  all  the  Lebesgue  integrals  are  well  defined  (this  hypothesis  holds 
true  for  p  €  %,  u,  v  G  V,  h  G  [L2(TAr)]d  1  and  f  G  [L2(fl)]d).  Then,  the  elastic 
problem  (3)  in  weak  form  reads,  for  a  given  material  distribution  p  G  "H: 

find  u  G  V  :  f?u( u;  p)(v)  =0  Vv  G  V,  p  G  H .  (5) 

We  observe  that  the  displacement  u  G  V  solving  Eqs.  (3)  and  (5)  depends  on  the 
prescribed  p,  for  which  u  =  u(p).  In  this  manner  the  stress  tensor  of  Eq.  (2)  associated 
to  the  solution  of  Eq.  (5)  reads: 

a(p)  :=  a(p ,  u(p))  =  C(p)e(u(p)).  (6) 

We  introduce  the  compliance  energy  of  the  system  (5),  say  Je{p),  as: 

Je(p)  ■=  f  V’e(p)  dfl,  (7) 

.In 

where  iI>e(p )  is  the  strain  energy  function: 

fpE(p)  '■=  i>E(p,u(p)),  (8) 

with: 

iPe{p,u)  :=  a(p,u)  :  e(u).  (9) 

We  observe  that  the  standard  definition  of  the  strain  energy  involves  a  factor  of  1/2 
which  is  neglected  to  maintain  the  formulation  consistent  with  the  typical  one  of  the 
topology  optimization  framework.  Also,  we  notice  that  following  from  Eqs.  (5)  and  (7), 
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JE(p)  =  /  f  •  u(p)  dfl+  ®  h  u(p)  dr  at-  Moreover,  Eq.  (7)  can  equivalently  be  writ- 

J  q  </  r  tv 

ten  as  Jb(p)  =  J£i(p,  u(p)),  where  from  Eq.  (9): 


j£;(p,u)  :=  /  'tjjE(p,u)dn. 

J  n 


(10) 


In  order  to  introduce  the  topology  optimization  problem  in  the  minimum  compliance 
case,  we  recall  that  only  a  limited  amount  of  material  can  be  used,  which  means  that 
only  a  limited  area/volume  of  Q,  say  V  <  |f2|,  can  be  covered  by  the  material,  with  V 
defined  as: 


V  :=  /  pd£l. 
Jn 


(ID 


In  this  sense,  the  topology  optimization  problem  is  constrained  and  the  space  of  ad¬ 
missible  controls,  say  Had  C  H,  in  which  we  look  for  the  optimal  solution  p*  is  defined 


as 


Had  :=  |</>  G  %  :  0  <  (j)  <  1  and  J  (j)  =  1.  Then,  the  problem  of  topology  op¬ 


timization  in  the  minimum  compliance  case  corresponds  to 


find  p*  G  Had  ■  P*  =  argmin  (  JE(p) )  (12) 

where  JE(p)  is  the  compliance  energy  (7)  of  the  elastic  system  (5). 

We  recall  that  the  elastic  properties  of  the  material  are  introduced  in  the  topology 
optimization  problem  by  means  of  the  interpolation  function  g(p)  of  Eq.  (1).  In  general, 
the  topology  optimization  procedure  allows  “soft”  regions  of  material  in  the  design 
domain,  since  p  can  assume  intermediate  values  between  0  (void)  and  1  (material). 
However,  these  situations,  even  if  consistent  with  the  mathematical  formulation,  are 
in  general  not  desired  as  output  of  the  optimization  problem.  Indeed,  the  goal  is  to 
distribute  a  given  material  with  its  full  elastic  properties  for  which  the  desired  values 
of  p  are  possibly  only  0  (void)  and  1  (material).  In  order  to  reduce  or  avoid  these 
situations,  the  function  g(p)  of  Eq.  (1)  plays  a  crucial  role.  In  the  material  interpolation 
formulation,  this  function  is  typically  chosen  such  that  intermediate  distributions  of 
material  correspond  to  a  material  with  poor  stiffness  properties;  in  this  manner,  due  to 
the  area/volume  constraint,  the  material  is  located  to  minimize  the  compliance  energy 
of  the  system.  A  typical  choice  for  g(p)  is  based  on  the  SIMP  model,  for  which: 

g(p)  =  pp,  (13) 


with  P  >  max  j - ,  — - l  if  d  =  2  or  P  >  max  j 15 - , - 1  if  d  =  3; 

\  1  -  i/0  1  +  i/0  /  l  7-5r0  21-2  v0  ) 

the  condition  on  the  power  P  is  to  guarantee  that  the  interpolation  model  represents  a 
material  model  [11,  12,  67].  Typically,  since  materials  with  i/o  =  1/3  are  often  consid¬ 
ered,  the  power  is  chosen  such  that  P  >  3  for  both  d  =  2,  3;  also  this  value  is  reasonable 
even  when  a  minimum  value  for  p,  say  pmin  >  0,  is  introduced.  We  observe  that  other 
choices  of  g(p)  can  be  made,  among  these  are  rational  functions  [99]  and  B-splines  [77], 
which  can  be  useful  for  vibration  problems;  see  [12]  for  a  wider  discussion. 

1  Tn  the  standard  definition  of  topology  optimization  problems  in  the  minimum  compliance  case,  the 
area/volume  constraint  is  an  inequality  one,  /  pdVl  <  V;  see  e.g.  [12].  However,  typically,  the  optimal 

Jn 

solution  p*  is  such  that  V *  :=  /  p*  dil  =  V  in  order  to  maximize  the  stiffness  of  the  structure. 

Jn 


The  continuous  topology  optimization  problem  (12)  does  not  admit,  in  general,  the 
existence  of  optimal  solutions  as  pointed  out  and  discussed  in  [12,  60];  moreover,  the 
uniqueness  of  the  solution  is  also  an  issue,  since  multiple  minima  can  be  detected  due  to 
the  non-convex  nature  of  the  problem.  However,  even  if  in  general  ill-posed  [60,  64,  95], 
the  topology  optimization  problem  is  discretized  and  then  solved  numerically.  Typically, 
see  e.g.  [11,  12,  99],  low  order  Finite  Element  approximations  are  used  to  approximate 
the  density  function  p,  even  if  other  choices  are  possible  (see  [45]  for  a  Finite  Volume 
Method).  Then,  the  problem  is  solved  by  means  of  a  suitable  optimization  method; 
in  this  sense,  the  Method  of  Moving  Asymptotes  (MMA)  represents  one  of  the  most 
effective  optimizers  for  the  solution  of  topology  optimization  problems  [100,  101].  The 
fact  that  the  continuous  problem  is  ill  -posed  reflects  on  the  numerical  solution,  even  if 
the  discrete  problem  admits  the  existence  of  optimal  solutions.  This  is  revealed  by  the 
so  called  mesh  dependency  effect,  for  which  different  optimal  solutions  are  obtained  with 
different  discretizations  of  the  problem,  specifically  for  different  Finite  Element  meshes. 
There  are  several  techniques  to  contain  or  eliminate  this  effect,  which  are  introduced 
as  global  or  local  constraints;  the  most  used  ones  are  [12,  95]:  local  constraints  on  the 
density  gradient  [79],  local  density  and  sensitivity  filters  [19,  79,  93,  94],  global  control 
of  the  minimum  length  scale  [49,  81],  perimeter  control  or  limitation,  for  which  a  global 
constraint: 


/  Vp 

Jn 


Vpdfl  <  Pl i 


(14) 


with  Pl  >  0,  is  imposed  [78].  The  imposition  of  a  perimeter  constraint  limits  the  number 
of  holes  that  a  solution  could  exhibit.  In  practice,  the  constraint  (14)  is  introduced  to 
make  the  topology  optimization  problem  well  posed,  as  shown  in  [7]  for  the  continuous 
case.  An  alternative  approach  to  the  inequality  constraint  on  the  perimeter  has  been 
proposed  in  [51]  and  consists  in  perturbing  the  objective  functional  (7)  with  a  smooth 
penalty  term  for  the  perimeter  Pl- 

Another  undesired,  but  common,  feature  that  optimal  topologies  can  exhibit  is  the  so 
called  checkerboard  phenomenon,  which  indicates  a  numerical  solution  with  patterns  of 
alternate  0-1  values  (void-material)  [12,  95].  As  discussed  in  [59],  this  represents  a  form 
of  numerical  instability  associated  with  the  approximation  of  the  topology  optimization 
problem,  which  is  a  nonlinear  mixed  variational  problem  in  the  independent  density  and 
displacement  variables.  Similarly  to  the  case  of  linear  mixed  problems  [22],  stability 
issues  could  arise  at  the  discrete  level  even  if  the  continuous  problem  is  well  posed. 
Even  if  a  comprehensive  analysis  of  the  stability  properties  has  not  been  performed  due 
to  the  nonlinear  nature  of  the  problem,  it  has  been  shown  in  [59]  for  some  particular 
topology  optimization  problems,  that  using  suitable  pairs  of  Finite  Element  spaces  for 
the  density  and  displacement  variables  eliminates  the  checkerboard  issues  2 .  In  general, 
as  pointed  out  in  [12],  any  of  the  numerical  techniques  introduced  to  limit  the  mesh 
dependency  effect  could  be  effectively  used  to  avoid  this  phenomenon. 

As  a  final  remark,  we  observe  that  in  order  to  properly  solve  the  topology  opti¬ 
mization  problem  (12)  in  the  SIMP  approach,  it  is  necessary  to  introduce  additional 
numerical  techniques,  which  in  general  depend  on  the  discretization  chosen,  with  re¬ 
spect  to  the  original  continuous  formulation  of  the  problem;  also  a  suitable  constrained 
optimizer  needs  to  be  used. 


2  Other  possibilities  are  the  introduction  of  an  augmented  Lagrangian  functional  or  the  postprocessing 
filtering  of  the  numerical  solution. 
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2.2  The  multiphase  approach 

In  order  to  introduce  the  topology  optimization  problem  in  the  multiphase  approach,  we 
observe  that  the  distribution  of  the  phases,  material  and  void,  inside  the  design  domain 
fl  is  represented  by  the  material  density  function  p.  Also,  we  define  the  parameters 
p  =  (7,  A)  £  V ,  with  7  and  A  positive  and  the  parameter  set  HcK2,  in  order  to  intro¬ 
duce  a  parametrization  on  the  topology  optimization  problem.  We  define  the  following 
penalized  objective  functional  which  depends  on  the  parameters  p  £  T>: 


J(pw)  ■=  /  i>{p'i  aO  dn, 

Jn 

(15) 

with  the  function  tp(p;  p)  defined  as: 

V>(p;  M)  ==  7 tMp)  +  Eo  (tMp)  +  A i>i(p)) , 

(16) 

where  the  strain  energy  function  iPe(p)  is  defined  in  Eq.  (8),  ij)B 
energy  function  and  ipi{p)  is  the  interface  energy  function: 

(p)  is  a  suitable  bulk 

II 

to  |  l-^ 
<1 

<1 

(17) 

The  constant  E0  is  introduced  to  ensure  that  all  the  terms  in  Eq.  (16)  have  the  dimension 
of  an  energy  density  3.  The  objective  functional  J(p;  p)  (15)  can  be  rewritten  as: 

J(p;  p)  =  JE{p]  7)  +  Jb(p\  Eq)  +  J/(p;  A;  E0) 

,  (18) 

where: 

Je(p;  7)  ■=  7  /  ipE{p)dn, 

Jn 

(19) 

Jb{p;E0 )  :=  E0  [  tpB{p)  dCl, 

J  n 

(20) 

Ji{p ;  A;  E0 )  :=  E0  A  f  tpi(p)  dfl. 

Jn 

(21) 

We  observe  that  the  objective  functional  (15)  is  penalized  in  the  sense  that  two  terms 
proportional  to  tpB(p)  and  ipi(p)  are  added  to  the  strain  energy  function  i/je(p)-  In 
particular,  the  bulk  energy  iI>b{p)  is  a  non-convex  smooth  function  chosen  in  the  form 
of  a  double-well  in  the  pure  phases  p  =  0  and  p  =  1,  for  example,  as  i/)b(p)  =  p(l  —  p)  or 
iPb(p)  =  p2{  1  ^  p)2  [23];  in  this  manner,  the  values  assumed  by  V’b(p)  for  intermediate 
values  of  p  are  larger  than  for  the  pure  phases,  which  are  preferred  in  the  optimization 
context.  Also,  further  penalization  terms  can  be  added  to  iPb(p)  in  proximity  of  the 
pure  phases  p  =  0  and  1  in  order  to  remove  the  inequality  constraints  0  <  p  <  1  from  the 
formulation  of  the  optimization  problem;  a  possibility  is  to  consider  iI>b{p)  a  logarithm- 
type  function  with  singularities  in  the  pure  phases.  The  interface  energy  function  ipi(p) 
plays  an  important  role  in  the  penalization  of  the  compliance  energy,  since  it  represents 
a  measure  of  the  perimeter  of  the  interfaces  between  the  phases;  in  this  sense,  it  is  the 
relaxed  version  of  the  global  perimeter  limitation  constraint  (14)  often  introduced  in 
the  SIMP  method  to  remove  the  mesh  dependency  effect  in  the  discretized  problem. 


^Eventually,  the  constant  Eq  could  be  included  among  the  parameters  fi  £  T>\  however,  the  para¬ 
metric  dependence  of  the  objective  functional  J(p;  //)  (15)  would  still  be  completely  represented  by  only 
two  parameters  by  means  of  suitable  scalings. 
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Also,  this  term  assumes  the  role  of  controlling  the  thickness  of  the  interfaces  through 
the  parameter  A  and  thus  the  capability  to  capture  the  geometrical  informations  from 
the  optimal  topology. 

Finally,  we  notice  that  the  objective  functional  (15)  can  be  rewritten  as  J(p;p)  = 
where: 

J \p,  u;  fi)  :=  [  ^{p,  u;p)dtt,  (22) 

Jn 

with,  following  from  Eq.  (9): 

${p,  u;  n)  :=  7 $e(p,  u)  +  E0  (ipB(p)  +  A tpi(p))  •  (23) 

Similarly  to  Eqs.  (18)— (21 )  we  can  write: 

J{Pi u;  n)  =  JE{p,  u;y)  +  JB(p;E0)  +J(p;\-E0),  (24) 

with: 

Je(p,  u;  7)  :=  7  /  i>E(p,u)dn.  (25) 

Jn 

If  we  consider  the  penalized  objective  functional  J(p;p)  (15)  with  the  bulk  energy 
function  V’b(p)  embedding  the  penalization  terms  for  the  inequality  constrains  0  <  p  < 
1,  we  can  take  the  space  of  admissible  controls  as  Had  =  {</>  £H  :  =  V},  where 

only  the  area/volume  constraint  is  explicitly  imposed  and  H  =  {<fi  £  n  U1(fl)  : 

g(cj))  e  L°°(n)};  the  extra  regularity  of  H  with  respect  to  the  SIMP  approach  is  due  to 
the  presence  of  the  interface  energy  ipi{p)  in  the  formulation.  The  topology  optimization 
problem  in  the  multiphase  approach  corresponds  to: 

find  p*  €  Had  ■  P*  =  argrnin  (  J(p;  /u) ) ,  (26) 

for  any  given  parameter  fi  £  T>.  We  observe  that  the  optimal  distribution  of  material  in 
the  design  domain  is  p*  =  p*  (x;  fi)  in  the  sense  that  it  depends  on  the  parametrization 
introduced  in  the  penalized  objective  functional  (15). 

The  parameters  p,  =  (7,  A)  £  T>  in  Eq.  (16)  play  a  crucial  role  in  the  optimization 
problem,  since  they  regulate  the  balance  of  the  different  terms  contributing  to  the  pe¬ 
nalized  objective  functional  J(p;/x).  For  example,  if  7  is  very  large,  the  optimization 
problem  assumes  a  similar  behavior  to  the  SIMP  approach.  Conversely,  if  7  =  0,  then 
the  optimal  result  p*  is  dominated  by  the  bulk  and  interfaces  terms  only;  this  last  case 
corresponds  to  a  pure  multiphase  problem,  where  the  total  free  energy  of  the  system  is 
minimized  [26,  40]  (see  Sec.  3  for  the  Cahn-Hilliard  equations). 

Remark  2.1  A  third  positive  parameter ,  say  k  £  R,  can  be  added  to  the  parameter 
vector  p,  £  V,  such  that  pK  :=  (7,  A,  n)  €  VK,  with  Cl3.  In  this  case,  the  penalized 
objective  functional  (15)  is  redefined  as: 


Up-,  ac) 


with  the  function  4>k,(p ;  Ac)-' 


iMp;aO  ■=  Ue{p)  +  e0 


+  kXipi(p) 


(27) 


(28) 
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Then,  the  formulation  of  the  topology  optimization  problem  follows  similarly  to  the  pre¬ 
vious  case.  The  role  of  k  is  twofold.  Firstly,  it  ensures  that  the  penalized  objective  func¬ 
tional  is  convex  for  k  “sufficiently”  large  and  the  topology  optimization  problem  (26)  is 
well-posed,  as  shown  in  [23]  for  a  minimum  weight  topology  optimization  problem  with 
stress  constraints  in  a  relaxed  approach.  Secondly,  it  ensures  that  the  optimal  topology 
converges  to  the  pure  phases  0  and  1  for  k  — >  0  [23],  according  to  the  properties  of 
r -convergence  (see  e.g.  [68])  for  functionals  with  interface  terms  [20].  On  this  basis  the 
parameter  k  introduces  the  possibility  to  solve  the  topology  optimization  problem  (26) 
by  means  of  the  continuation  method  [12,  20,  23,  71].  In  this  procedure,  the  optimal 
topology  p*  is  obtained  as  the  last  step  of  a  sequence  of  minimizers  of  locally  convex 
(or  quasi-convex)  optimization  problems  parametrized  for  decreasing  values  of  k  and 
initialized  with  the  optimal  solution  of  the  previous  step;  this  means  that  a  non-convex 
topology  optimization  problem  without  interface  term,  which  corresponds  to  the  ideal 
formulation,  is  obtained  as  the  limit  for  k  — >  0. 

Finally,  we  observe  that  the  topology  optimization  problem  in  the  multiphase  ap¬ 
proach  is  completely  defined  at  the  continuous  level  and,  eventually,  for  suitable  choices 
of  the  parameters,  also  well  -posed.  This  is  not  the  case  of  the  SIMP  formulation,  which 
is  in  general  ill-posed  at  the  continuous  level  and  it  requires  the  introduction  of  addi¬ 
tional  numerical  techniques  at  the  discrete  level.  However,  a  non-convex  optimizer  is 
still  required  to  solve  the  topology  optimization  problem  in  the  multiphase  context. 


3  Phase  Field  Model:  the  Cahn— Hilliard  Equation 

In  Sec.  2.2  we  have  considered  the  topology  optimization  problem  as  a  multiphase  ap¬ 
proach  for  which  a  penalized  objective  functional  is  minimized;  the  formulation  is  how¬ 
ever  still  set  in  an  optimal  control  context,  for  which  an  optimization  problem  needs 
to  be  solved  to  find  p* .  However,  the  penalized  objective  functional  can  be  interpreted 
as  a  total  free  energy  and  the  topology  optimization  problem  recast  in  a  phase  transi¬ 
tion  setting.  Indeed,  in  this  case,  the  evolution  of  the  phases  is  such  that  an  energy  is 
minimized  in  time  with  respect  to  the  initial  configuration.  In  this  section  we  provide 
the  derivation  of  the  standard  Cahn-  Hilliard  phase  field  model  [25,  26,  27,  28]  and  we 
highlight  its  properties  in  view  of  the  phase  field  model  for  topology  optimization  in 
Sec.  4.  For  further  details  in  the  analysis  of  the  Cahn-Hilliard  equation  we  refer  the 
interested  reader  to  [24,  37,  38,  40,  85],  while  e.g.  to  [37,  46,  103]  for  its  numerical 
solution. 

Let  us  consider  a  two  phase  problem,  for  which  the  phase  transition  is  described  by 
the  variable  p  =  p(t,x),  which  corresponds  to  the  concentration  of  one  of  the  phases 
in  the  domain  Q  (the  other  one  is  obtained  as  1  —  p  in  the  0-1  binary  representation). 
The  total  free  energy  (Ginzburg-Landau  free  energy),  which  we  indicate  with  F(p-,  A) 
[26,  27,  40],  is  expressed  in  the  case  of  interest  as: 

F(p-,X):=  f  ifF{p;\)dn,  (29) 

Jn 

with  the  total  free  energy  function  i]>F(p ;  A)  defined  as: 

ipF(p;  A)  :=  C0  {ipB(p)  +  M Pi(p))  ,  (30) 

for  a  positive  parameter  AsPcR  and  the  bulk  and  interface  energies  given  in  Sec.  2.2; 
the  parameter  Co  is  introduced  to  ensures  that  the  function  ipF(p;  X)  assumes  the  di¬ 
mension  of  an  energy  density.  Typically,  a  double-well  logarithm  or  quartic  function  is 
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chosen  for  V’b(p)  [26,  38,  40,  103],  even  if  other  type  of  functions  can  be  conveniently 
used. 

The  definition  of  the  phase  transition  model  is  based  on  the  concept  of  gradient  flow 
gradF(p,  A)  of  the  functional  F(p;  A)  in  the  norm  of  a  Hilbert  space  Z  [28,  40],  for  which 
the  corresponding  equation  in  strong  form  reads: 

^  =  -grad  F(p;  A)  in  ft,  Vt  G  [0,  T),  (31) 

with  p  =  po  for  t  =  0  in  ft.  Depending  on  the  choice  of  the  function  space  Z ,  different 
phase  transition  models  can  be  obtained;  if  Z  =  j<p  G  (i^1(ft)),  :  (ip,  1)  =  0  j  [40],  the 

Cahn-Hilliard  gradient  flow  and  the  corresponding  equation  are  obtained  [25,  26,  27], 
otherwise,  if  Z  =  L2(f2),  the  Cahn-Allen  equation  is  derived  [6].  In  particular,  in  the 
Cahn-Hilliard  case,  the  gradient  flow  grad  F(p,  A)  reads: 

grad  F(p;  A)  =  -V  •  ( M(p)WzF(p ;  A)) ,  (32) 

where  M (p)  >  0  is  a  sufficiently  regular  function  called  the  mobility,  which  is  typically  a 
constant  M0  or  degenerate  M(p)  =  M0p(l— p),  and  zF(p;  A)  is  the  potential  associated  to 
the  total  free  energy  function  ipF(p ;  A);  in  the  Cahn-Allen  case,  the  gradient  flow  would 
read  gradF(p;A)  =  M(p)zf(p',  A).  The  potential  Zf(p\  A),  which  we  introduce  for  the 
sake  of  simplicity,  depends  specifically  on  the  choice  made  for  the  functional  (29)  and 
in  this  case  is  obtained  as  the  Gateaux  derivative  in  L2(ft)  of  the  total  free  energy  (29): 

1  dF 

(zF(p-,\),(j>)L2{u)  =  —  —  {p-,\)[<j>\  V^GU'fO),  (33) 

for  some  A  G  T>\  if  further  we  assume  that  the  boundary  condition  Vp  •  n  =  0  is  imposed 
on  clft,  the  potential  Zp(p]  A)  simply  reads: 

zF(p;X)  =  zB(p)  +  Xzi(p),  (34) 


where: 

,  ,  dtps,  , 

zb(p)  ■= 

dp 

Zi(p )  :=  —A p, 

for  p  sufficiently  regular. 

A  standard  formulation  of  the  Calm  Hilliard  equation  in  strong  form  is: 
§£  =  V  •  (. M(p)\7zF(p ;  A))  in  ft,  Vt  G  [0,  T), 


M(p)VzF(p\  A)  •  n  =  0 
Vp  •  n  =  0 
P  =  Po, 


on  dfl,  Vt  G  [0,T), 
on  dfl,  Vt  G  [0,T), 
in  ft,  t  =  0, 


(35) 

(36) 


(37) 


for  some  A  G  V.  We  define  the  function  space  T~L  :=  {(p  G  if2(fl)  :  V</>  •  n  =  0}.  Let  us 

assume  that  p  G  C1  ([0,  T);  H)  and  so  ^  G  C°  ([0,  T);  H).  That  is,  p  is  a  C'1-continuous 

dp 

mapping  from  the  time  interval  [0,  T)  into  TL  and  '  is  a  C° -continuous  mapping  from 
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[0 ,  T)  into  H.  Consequently,  for  each  t  G  [0 ,  T),  p  e  H  and  —  G  H.  From  Eq.  (37),  for 
each  t  G  [0 ,  T),  the  residual  Rch{p ;  G  M  is  given  by: 


Rch{P'A){4>)  :=  [  ^<j>dSl+  l  M(p)S7zF(p;\) -V^dCl.  (38) 

Jn  ot  Jn 

Then,  the  Cahn-Hilliard  equation  in  weak  form  reads: 

find  :  RCH(p;\){<1>)  =  0  V</>  €  "H,  Vf  <E  [0,  T), 

(39) 


with  p  =  po  in  fi,  t  =  0, 

for  some  A  €  V.  The  Cahn-Hilliard  equation  (39)  (or  Eq.  (37))  is  endowed  with  the 
following  properties: 


•  Its  solution  p  exists  and  is  unique  for  the  case  of  constant  mobility  as  shown  in 
[85]  for  problems  in  dimensions  d  =  1,2,3  under  suitable  hypothesis  for  the  bulk 
function  i/)b(p)  including  its  smoothness;  existence  of  solutions  is  discussed  in  [38] 
for  the  case  of  degenerate  mobility. 

•  It  is  mass  conservative  in  the  sense  that  the  area/volume  covered  by  the  phases  in 
ft  is  constant  in  time;  indeed,  by  using  the  definition  (11),  we  can  easily  deduce 
from  Eq.  (39)  for  <j)  =  1  that: 


Vi  €  [0,  T). 


(40) 


•  The  total  free  energy  functional  (29)  is  a  Liapunov  functional;  indeed,  it  is  possible 
to  show  from  Eq.  (39)  that: 


dF  f 

—  (p;  A)  =  -Co  /  M(p)V2F(p;  A)  •  VzF(p-:  A)  dft  <  0  Vi  G  [0,  T).  (41) 

at  J  a 

This  implies  that  the  phase  transition  occurs  in  such  a  manner  that  the  energy 
associated  to  the  Cahn-Hilliard  equation  is  decreasing  or  at  most  conserved  in 
time;  this  property  also  holds  true  for  the  Calm  Allen  equation,  since  in  general: 


—  (p;A)  =  — Co||grad F(p;  A)|||  <  0, 


(42) 


in  the  norm  induced  by  Z ,  see  [40]. 

•  Under  suitable  hypothesis  on  the  function  i/jb{p)  including  its  analyticity,  in  [85] 
it  is  proved  that,  for  a  given  po,  the  unique  solution  p  converges  to  an  equilibrium 
dp 

(steady  state)  and  —  — >  0  for  t  — >  oo  in  the  topology  of  the  corresponding  function 

spaces.  This  implies  from  Eqs.  (31)  and  (42)  that  the  steady  state  solution  is  a 

dF 

critical  one  for  the  total  free  energy  F(p;  A)  [37],  with  — (p;  A)  — >  0  for  t  — >  oo;  it 

follows  that  F(p\  A)  evolves  to  a  local  minimum  through  the  phase  transition  from 
the  initial  solution  po- 
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4  Topology  Optimization  with  the  Phase  Field  Model 

We  derive  now  the  phase  field  model  for  topology  optimization  similarly  to  [114].  First, 
we  provide  the  generalized  Calm -Hilliard  equation  based  on  the  multiphase  approach 
of  Sec.  2.2;  then,  we  reformulate  the  problem  as  a  coupled  system  with  phase  and  dis¬ 
placements  as  independent  variables  and,  finally,  we  discuss  the  dimensionless  problem 
highlighting  its  dependence  on  dimensionless  parameters. 


4.1  The  generalized  Cahn-Hilliard  equation 

In  Sec.  3  we  derived  the  Cahn-Hilliard  equation  starting  from  a  total  free  energy  func¬ 
tional.  We  have  observed  that  this  phase  field  model  is  area/volume  conservative  and 
the  phase  transition  occurs  in  such  a  manner  that  the  energy  of  the  system  decreases  in 
time.  These  two  properties  are  crucial  to  recast  the  multiphase  approach  for  topology 
optimization  of  Sec.  2.2  in  a  phase  transition  model,  since  an  area/volume  constraint  is 
set  in  the  minimum  compliance  case  and  a  penalized  objective  functional  needs  to  be 
minimized. 

The  derivation  of  the  generalized  Cahn-Hilliard  equation  for  the  topology  optimiza¬ 
tion  problem  in  the  minimum  compliance  case  follows  in  similar  manner  to  Sec.  3  by 
using  the  penalized  objective  functional  (15).  In  particular,  we  have  that: 

dp 

—  =  -grad  J(p;n)  in  O,  Vi  €  [0,  T),  (43) 

where  p  =  po  for  t  =  0  in  Q, 

grad  J(p;/d  =  -V  •  (M (p)V z(p\  p,)) ,  (44) 


and  the  potential  z(p\  p)  deduced  from  the  Gateaux  derivative  in  L2(fl)  of  the  penalized 
objective  functional  (15): 


Wp;m).  4>)l*(si)  =  v)W\  v^etf1^),  (45) 

for  some  /i€D.  If  we  assume  that  Vp  •  n  =  0  on  <9f 2,  we  have  from  Eq.  (16)  that: 

z(p;p,)  =  ^rzE{p)  +  zB(p)  +  XzI(p),  (46) 

where  Zb(p)  and  Zi(p)  are  defined  in  Eqs.  (35)  and  (36),  respectively,  and: 

Ze(p)  ■=  d^{p)-  (47) 

dp 


In  order  to  evaluate  Ze(p),  further  elaborations  are  needed  since 
from  Eq.  (8);  in  particular,  we  have  that: 


dips 

dp 


(fou(p)) 


ze(p) 


dipE 

dp 


(p,u(p))  + 


dipE 

d\i 


(p,u(p)) 


(48) 


From  Eq.  (9)  we  deduce  that: 


dipE 

dp 


(P, u) 


dd 

dp 


(p,u)  :  e(u); 


(49) 
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similarly,  by  recalling  that  a(p,  u)  depends  linearly  on  u  and  the  elastic  tensor  C (p)  is 
symmetric,  we  have: 


dips 

du 


(fou) 


du 

dp 

P,  t:  )  :  e(u)  +cr(p,u)  :  £  (  ^ 


=  2a  p, 


du 


du 

dp 

du 

dp 


(50) 


:  £(u). 


In  order  to  evaluate  the  term  a  p.  —  ,  we  need  to  differentiate  the  weak  form  of  the 

V  dpj 

elasticity  equation  (5)  with  respect  to  p;  by  assuming  that  the  function  —  G  V,  we 

dp 

obtain:  __ 

J  d  (p,  :  e(v)dfl  +  J  ^(p,  u)  :  e(v)  dfl  =  0  Vv  G  V,  (51) 

and  hence  for  v  =  u: 


a 


£(u) 


dd 

dp 


(p,u)  :  £(u) 


dips 

dP 


(fou). 


(52) 


By  replacing  the  result  (52)  in  Eq.  (50),  and  then  Eqs.  (49)  and  (50)  in  Eq.  (48),  we 
obtain: 

zE(p)  =  zE(p,u(p)),  (53) 

with:  __ 

2e(p,u)  :=  -^^(p,u).  (54) 

The  potential  (46)  can  also  be  written  as: 

z(p\n)  ■=  ^(P,u(p);pt),  (55) 


where: 


and  equivalently: 


z(p,  u;  /x)  :=  ~^zE{p,  u)  +  zB(p)  +  A z7(p), 


x  7  dxpE ,  ,  ,  dips,  ,  .  A 

z[P,  u;m)  :=  -vrqr-(P,u)  +  -7—  (p)  -  AAp. 


E0  <9p 


dp 


(56) 


(57) 


It  is  now  possible  to  introduce  the  strong  form  of  the  generalized  Cahn-Hilliard 
equation  for  topology  optimization  from  Eq.  (43),  which  reads: 


^=V-(M(p)V*(p;/*)) 

M(p)S7z(p-,n)  •  n  =  0 
Vp  •  n  =  0 
P  =  Po  > 


in  II,  Vi  G  [0, T), 
on  dCl,  Vi  G  [0,  T), 
on  c?0,  Vi  G  [0,  T), 
in  II,  i  =  0, 


(58) 


for  some  p,  £  V.  By  recalling  from  Sec.  3  that  p  G  W  for  all  t  G  [0,T),  with 
H:=  {(p£  H2(fl)  :  V0  •  n  =  0 } ,  we  introduce  the  residual  Rp(p',  G  R.  such  that: 


RP{p\P »)(</>)  ■  =  J  ^  (pdn  +  M(p)Vz(p;n)  -S7(pdCl. 


(59) 
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It  follows  that  the  weak  form  of  the  generalized  Cahn-Hilliard  equation  is: 

find  p  €  W  :  Rp(p(t;  p);  p) (</>)=  0  V</>  €  H,  Vf  G  [0,  T), 

with  p  =  po  in  12,  t  =  0, 


(60) 


for  some  p  £  V.  We  observe  that  Eq.  (60)  is  obtained  with  the  natural  boundary 
condition  M(p)Vz(p ;  /x)  -  n  =  0  and  the  essential  one  Vp-  n  =  0  defined  on  the  boundary 
(90;  due  to  its  nature,  the  latter  is  embedded  in  the  space  T~L. 

The  generalized  Cahn-Hilliard  equation  (60)  (or  Eq.(58))  represents  a  model  for 
topology  optimization  problems  in  the  minimum  compliance  case,  for  which  the  following 
properties  hold  similarly  to  the  Cahn-Hilliard  equation: 

•  The  unique  solution  p  exists  by  extending  the  result  of  [85]  under  suitable  hy¬ 
pothesis  on  the  bulk  and  strain  energy  functions  iPb(p)  and  ^e(p),  in  the  case  of 
constant  mobility. 

•  The  area/volume  covered  by  the  material  in  the  design  domain  12  is  constant 
during  the  phase  transition;  see  Eq.  (40). 

•  The  penalized  objective  functional  J(p;p)  (15)  is  a  Liapunov  functional  which 
evolves  in  time  by  decreasing  from  the  initial  value  corresponding  to  po;  indeed, 
from  Eq.  (42)  in  analogy  with  Eq.  (41),  we  have: 

7  t  n 

=  ~EoJ^  M(p)Vz(p;p)  ■  Vz(p-,p)dfl  <  0  \/fe[0,T).  (61) 

•  If  the  functions  i/)b(p)  and  ?/2e(p)  satisfy  the  hypothesis  made  in  [85]  for  the  bulk 
energy  of  the  Cahn-Hilliard  equation,  the  unique  solution  p  converges  to  a  steady 
state  which  minimizes  the  penalized  objective  functional  J(p;  p)  of  the  multiphase 
approach  with  respect  to  the  initial  solution  po;  in  this  manner,  the  optimization 
problem  is  converted  to  a  time  dependent  one  and  its  optimal  solution  p*  is  ob¬ 
tained  as  p  for  t  — >  oo. 

•  The  topology  optimization  problem  is  completely  defined  in  the  formulation  by 
choosing  the  data  h,  f,  Co  and  12,  the  function  g{p)  in  Eq.(l),  the  mobility  M(p), 
the  initial  condition  po  (which  also  introduces  the  area/volume  constraint  V )  and 
the  parameters  p  £  T>. 

In  the  current  work  we  do  not  provide  a  rigorous  analysis  of  the  generalized  Cahn- 
Hilliard  equation  for  topology  optimization  problems.  At  this  point,  we  only  speculate 
on  the  possibility  that  the  bulk  iPb(p)  and  especially  the  strain  energy  iPe(p)  functions 
could  satisfy  the  hypothesis  made  in  [85]  for  the  existence  and  uniqueness  of  the  solution 
p  and  its  convergence  to  a  steady  state,  the  minimizer  of  J(p'1p);  further  limitations 
on  the  material  interpolation  model  could  be  eventually  deduced  to  fit  such  hypothesis. 
However,  we  observe  that  numerical  tests  exhibit  the  convergence  of  the  solution  to  a 
steady  state  for  any  given  initial  solution,  for  which  we  feel  that  the  validity  of  the 
considered  phase  field  model  could  be  shown  also  from  an  analytical  point  of  view. 


4.2  The  coupled  system 

The  generalized  Cahn-Hilliard  equation  (60)  (or  Ecp  (58)  in  strong  form)  is  written  in 
terms  of  the  phase  variable  p.  However,  in  order  to  solve  the  problem  it  is  necessary  to 
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evaluate  the  displacement  u(p)  by  solving  the  elasticity  equation  (5)  for  a  given  value 
of  p.  In  practice,  for  each  t  €  [0,T),  there  is  a  corresponding  displacement  u  €  V. 
For  this  reason,  it  is  convenient  to  consider  the  phase  p  and  the  displacement  u  as  two 
independent  variables  and  rewrite  the  generalized  Cahn-Hilliarcl  equation  as  a  coupled 
system  of  equations. 

By  recalling  Eqs.  (58),  (56)  and  (3),  the  strong  form  of  the  coupled  system  is: 


^  =  V  •  (M(p)W(p,  u;  /*)) 
-V  •  d{p,  u)  =  f 
M(p)Vz{p,  u;n)  •  n  =  0 
Vp  •  n  =  0 

u  =  0 

a(p,  u)n  =  h 
P  =  Po, 


in  Cl,  Vi  G  [0,  T), 
in  SI,  Vi  G  [0,  T), 
on  DCl,  Vi  G  [0,  T), 
on  dCl,  Vi  G  [0,  T), 
on  rD,  Vi  G  [0,  T), 
on  VN,  Vi  G  [0, T), 
in  Cl,  t  =  0, 


(62) 


for  some  fi  G  V.  Also,  for  each  t  G  [0,  T),  we  define  the  residual  Rp(p,  u;/Lt)(0)  G  R 
from  Eq.  (59)  by  recalling  the  potential  z(p,u;p,)  (56): 

RP{p,^',  P){4>)  :=  j  ^  <j>dCl  +  M(p)Vz(p,u;n)  -C/tpdCl-  (63) 

similarly,  from  Eq.  (4),  for  each  t  G  [0,  T),  we  redefine  Ru(p,  u;  /x)(v)  G  R  to  highlight 
the  explicit  dependency  on  p  and  u  as: 


Ru(p,  u;  |i)(v)  :=  j  <j{p,  u)  :  e(v)  dCl  —  f  f-vdfl—  ®  h-vdrAr.  (64) 
J O  J  £1  v  Tjv 

Then,  the  coupled  system  of  generalized  Calm  -Hilliard  and  elasticity  equations  in  weak 
form,  which  we  will  indicate  as  TO(/i),  reads: 

find  p  G  W,  u  G  V  : 

Rp(p,  u;*i)(0)  =  O  V^GH,  ViG  [0,T) 

TO(/x)  _  (65) 


Ru(p,  u;  /Lt)(v)  =  0  Vv  G  V,  Vi  G  [0,T), 

with  p  =  po  in  fl,  i  =  0, 


for  some  fi  G  2?  and  the  associated  energy  J(p,  u;/x)  defined  in  Eq.  (22).  We  observe 
that  the  displacement  u  depends  on  the  time  f  G  [0,  T)  implicitly  trough  the  variation 
of  the  phase  p.  Indeed,  in  the  elasticity  equation  no  time  derivatives  appear,  meaning 
that  the  displacement  adapts  instantaneously  to  the  variation  of  the  phase. 

The  reformulation  of  the  generalized  Cahn-Hilliard  equation  (60)  (or  Eq.  (58))  into 
the  coupled  system  TO(/x)  (65)  provides  an  approach  to  analyze  the  problem  in  terms  of 
existence  and  uniqueness.  Indeed,  we  notice  that  the  coupled  system  (65)  shows  many 
analogies  with  the  Cahn-  Hilliard  equation  with  elastic  misfit,  also  known  as  the  Cahn- 
Larche  model  [62],  for  which  the  phase  transition  is  also  driven  by  elastic  interactions 
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of  the  material  [42,  73].  In  [42]  the  existence  of  a  solution  of  such  a  system  is  proved  as 
well  as  its  uniqueness  for  a  specific  choice  of  the  elastic  energy;  in  [43]  the  corresponding 
discretized  problem  is  analyzed.  However,  we  remark  that  the  Cahn-Larche  equations 
and  the  coupled  system  TO(/i)  (65)  corresponding  to  the  generalized  Cahn-Hilliard 
equations  are  derived  from  different  concepts.  Indeed,  the  first  equations  follow  from 
thermodynamical  considerations  and  balance  laws  for  the  species  and  the  momentum 
with  p  and  u  as  independent  variables.  Conversely,  the  generalized  Cahn-Hilliard  equa¬ 
tion  is  developed  by  considering  the  displacement  variable  as  dependent  on  the  phase 
variable  u(p)  and  only  the  balance  of  the  species  is  taken  into  account;  the  coupled 
system  TO(/i)  only  represents  a  reformulation  of  such  an  equation. 


4.3  Dimensionless  form  of  the  coupled  system 

We  now  rewrite  the  coupled  system  (65)  in  dimensionless  form.  With  this  aim,  we 
introduce  the  dimensionless  space  and  time  coordinates: 

x*  =  Xj/ Lq,  i  =  l,...,d,  t*  =  t/T0  (66) 

and  the  phase  and  displacement  variables: 

P *  =  P,  u*  =  u/L0,  (67) 


where  Lq  and  To  are  representative  length  and  time  scales,  while  the  superscript  *  indi¬ 
cates  dimensionless  variables.  Also,  if  we  use  Eq  as  the  representative  Young  modulus, 
we  obtain: 

d*(p*,u*)  =  a(p,  u )/E0, 


e*(u*)  =  e(u), 


and 


z*E(p*,  u*)  =  zE{p ,  u)/ To, 

=  zB(p), 

4(4)  =  Lo  zi(p )> 

while  defining  the  reference  surface  and  body  forces  Hq  and  /o,  we  have: 

h*  =  h//So,  f*  =  f//0; 


u*)  =i>E{p,v)/EQ, 

Vb(P*)  =  i>B(p), 

VAp *)  =  Lli>i{p), 


(68) 


(69) 


(70) 


finally,  the  dimensionless  mobility  is: 


M*(p*)  =  M(p)/M0. 


(71) 


Let  us  define  the  following  dimensionless  parameters  which  we  indicate  with  the  vector 
D  =  (Di, . . . ,  D5)  such  that: 


D\  := 
D4  := 


T  4 
^0 


T0AM0’ 

Eo 

h’ 


D  -Ll 

-  T, 


03:= 4 


Da  := 


foLo 

^T; 


(72) 


we  observe  that  the  parameter  7  is  dimensionless,  while  the  parameter  A  assumes  the 
same  dimension  as  Lq.  Moreover,  we  define  the  dimensionless  potential: 


2* (4,  u*;  D)  :=  D3  z*E(p\  u*)  +  D2  z*B{p*)  +  z*j(p*). 


(73) 
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for  which  z*{p *,  u*;  D)  =  D2z(p ,  u;  fi)  and  the  dimensionless  energy  function: 

r(p\  u*;  D)  :=  D3  frE{p* ,  u  *)  +  D2  V£(P*)  +  #(P*)»  (74) 

D2  ~ 

for  which  ip*(p*  ,u*;D)  =  — ip(p,u;  p,).  With  these,  we  define  from  Eq.  (63)  the  di- 
__  E0 

mensionless  residual  u*;D)(0*)  eMVle  [0,  T): 


i?7(p*,u*;D)(0*)  :=  Dx  I  ~-£-  <p*  dfi* 


dp* 

at* 


+  [  M*(p)V*?*(p*,  u*;  D2j  D3)  •  V*</>*  dfi* 
Jn* 


(75) 


o 

4— d 
0 


for  which  R*{p*,  u*;  D)((/>*)  =  u  Rp(p,u;  n){<j>)\  similarly,  from  Eq.  (64)  we  define 


-Pyr  /  -  AMq  -P 


I?*(p*,u*;D)(v*)  G  K.  Vt  G  [0,  T)  as: 

i?*(p*,u*;D)(v*)  :=  D4  [  5*(p*,  u*)  :  £*(v*)  dfi* 


-£i5  /  f*  •  v*  dfi*  -(f)  h*  •  v*  dT^, 


(76) 


with  i?*  (p*,  u*;  D)(u*)  =  7 — 7  i?u(p,  u; /x)(v).  It  follows  that  the  dimensionless  topo- 

^0-^0 

logy  optimization  problem  in  the  phase  field  approach  reads: 

find  p*  G  W,  u*  G  V  : 
i?*(p*,u*;D)(<A*)=0 
I?*(d*,u*;D)(v*)  =  0 
with  p*  =  Pq 


TO*  (D) 


W  G  R,  Vt*  G  [0,  T*) 
Vv*  G  V,  Vt*  G  [0,  T*), 
in  fi*,  t*  =  0, 


(77) 


with  p^  =  pq  and  T*  =  T/Tq.  The  corresponding  dimensionless  penalized  objective 
functional  (energy)  follows  from  Eq.  (74)  as: 


J*(p*,u*;D)  =  [  i/j*(p*,  u*;  D2,D3)  dfi*, 
Jn* 


(78) 


~  1  L2~d  ~  ~ 

with  J*(p*,  u*;  D)  =  — - ^ — J(p,  u;  /x);  the  same  scaling  occurs  for  J%(p* ,  u*),  J%{p*) 

hjQ  A 


and  Jj(p*). 

We  notice  that  the  topology  optimization  problem  (77)  fully  depends  on  the  dimen¬ 
sionless  parameters  D  and  the  initial  solution  p^,  for  which  the  principle  of  dimensional 
similitude  can  be  applied.  Also,  it  is  possible  to  reduce  the  parametric  dependence  by 

L4 

setting  Di  =  1  and  hence  choosing  the  representative  time  scale  as  T0  =  — — — ;  it  follows 

AMq 


that  the  problem  TO*  (D)  (77)  depends  on  D  =  (D2, . . . ,  D5). 

We  observe  from  Eq.  (76)  that  the  dimensionless  parameters  D4  and  D5  of  Eq.  (72) 
only  affect  the  linear  elasticity  equation.  The  dimensionless  parameters  D2  and  D3  are 
responsible  for  changing  the  optimal  solution  of  TO*  (D)  once  the  data  are  set;  indeed 
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they  control  the  balance  of  the  strain,  bulk  and  interfaces  energies  in  the  penalized 
objective  functional  (78).  From  Eq.  (74),  we  deduce  that  if  D3  is  very  large,  then  the 
phase  field  solution  is  dominated  by  the  minimum  compliance  term  of  the  penalized 
objective  functional.  Conversely,  if  D3  — >  0,  specifically  if  due  to  7  — >  0,  TO*  (D)  solves 
the  standard  Cahn-Hilliard  equation  for  the  phase  problem  with  the  displacement  u* 
depending  on  the  distribution  of  p*  in  Cl*.  Also,  since  from  Eq.  (72)  D3  =  7D2  and 
D-2  ~  1/A,  with  A  related  to  the  control  of  the  thickness  of  the  interfaces,  we  have  that 
D-2  and  D3  1,  when  sharp  interfaces  are  required. 

Remark  4.1  For  the  sake  of  simplicity,  we  henceforth  omit  the  superscript  *  to  indicate 
dimensionless  quantities. 

5  Numerical  Approximation 

In  this  section  we  discuss  the  numerical  approximation  of  the  topology  optimization 
problem  TO(D)  (77)  by  using  a  similar  approach  to  the  one  used  in  [46,  48]  for  phase 
field  problems.  In  particular,  for  the  spatial  approximation  we  consider  Isogeometric 
Analysis  [33,  55],  while  for  the  temporal  approximation  we  use  the  generalized-a  method 
[31]  with  time  adaptivity. 

5.1  The  spatial  approximation 

We  consider  design  domains  Cl  described  by  a  NURBS  (or  B-spline  basis)  [80],  for  both 
the  two  and  three-dimensional  cases.  The  use  of  Isogeometric  Analysis  for  the  spatial 
approximation  of  the  PDEs  allows  us  to  encapsulate  directly  the  geometry  representation 
in  the  analysis,  by  using  the  same  basis  functions  used  to  represent  the  geometry  [33,  55]. 
In  this  manner,  no  geometrical  approximations  are  introduced  in  the  analysis  of  the 
topology  optimization  problem  TO(D)  (77). 

Isogeometric  Analysis  provides  a  way  to  easily  achieve  high-order  continuity  in  the 
approximated  solution  without  introducing  extra  degrees  of  freedom.  In  particular,  for 
the  TO(D)  problem,  the  use  of  globally  C1  (O)-continuous  basis  functions  is  necessary 
to  approximate  the  functional  space  77  of  the  phase  variable  p.  Specifically,  we  consider 
basis  functions  of  degree  p  >  2  defined  by  open  knot  vectors  with  equally  spaced  internal 
knots  repeated  at  most  p  —  1  times;  in  this  manner  we  ensure  that  the  basis  functions 
are  globally  C,1(fl)-continuous.  We  consider  the  same  basis  functions  for  the  phase  p 
and  the  displacement  u. 

We  introduce  the  finite  dimensional  spaces  77/,  C  77  and  V/,  C  V,  with  the  above 
mentioned  properties,  of  dimensions  n/,iP  :=  dim  (77/,)  and  n/,)U  '■=  dim  (V/,).  Due  to  the 
properties  of  the  NURBS  basis,  the  essential  boundary  condition  Vp  ■  n  =  0  is  easily 
introduced  in  the  space  77/,  by  imposing  equality  of  two  consecutive  control  values  of 
p  normal  to  the  boundary.  We  illustrate  strong  imposition  of  the  essential  boundary 
condition  in  Fig.  2  for  a  one-dimensional  B-spline  basis  of  degree  2.  This  also  holds  for 
NURBS  and  for  the  multidimensional  case.  For  each  t  £  [0,T),  p/,  £  Hh  and  u /,  £  V/, 
have  the  representations: 

nbf  nbf 

Ph(t,  x)  =  y;  pA(t)NA(x),  u/,(f,  x)  =  y  u A(t)NA(x),  (79) 

A=1  A=1 
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\ 

Figure  2:  B-spline  basis  of  degree  2  and  knot  vector  2  = 

{{0}f=1, 0.2, 0.4, 0.6, 0.8,  {l}f=1};  the  two  external  pairs  of  basis  functions  are 
marked  with  squares  and  circles  to  indicate  that  the  corresponding  control  variables 
need  to  be  equal  to  impose  the  boundary  condition  Vp  ■  n  =  0  on  dfl. 


with  NA(x)  the  NURBS  basis  and  nbf  the  number  of  basis  functions;  the  corresponding 
test  functions  are: 

nbj  nbf 

<Ah(x)  =  <j)ANA(x),  vh(x)  =  ^2  vANA{x).  (80) 

A= 1  A=1 

Then  the  discrete  TO^(D)  coupled  problem  for  any  t  £  [0,T)  is: 

find  ph  £  Hh ,  e  Vh  ■ 

TOh(D)  R„(ph,  ufc;D)(^h)=0  Vfo  G  Uh,  (81) 

^u(Ph,Uh;D)(vh)  =  0  Vvfc  £  vh, 

where  the  discrete  initial  condition  po,/i  is  obtained  as  the  L2(Q)  projection  of  po  in  the 
space  Hh-  The  total  number  of  degrees  of  freedom  of  TO^(D)  is  nh  :=  nh,P  +  nh,u- 
Finally,  we  introduce  from  Eq.  (79)  the  vectors  of  control  variables: 

P(t):={pA(t)}^lv  P(t)  :=  {W(C,  U(t)  :=  {uA(t)}^x ,  (82) 

and,  from  Eqs.  (75)  and  (76),  the  discrete  residuals: 

Rp  (p(t),P(i),U(i))  :=  {i?p(pft,u/l;D)(lVA)}'^i, 

{r  ~  id  t  "*>/  (83) 

{i^/,,ufc;D)(JVAej)}.=iJ  , 

where  e$,  with  i  =  1, . . .  ,d,  represent  the  orthonormal  basis  of  the  space  for  the 
sake  of  simplicity,  we  neglected  the  explicit  dependency  of  the  discrete  residuals  on  D 
in  Eq.  (83). 
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5.2  The  time  approximation 

The  time  approximation  of  the  TO(D)  problem  (77)  represents  a  challenge  similar 
to  that  for  the  Cahn-Hilliard  equation  due  to  the  fourth-order  term  and  significant 
nonlinearities.  We  use  the  generalized-a  method;  see  [31,  57]  and  also  [8,  46].  Moreover, 
since  the  solution  of  the  topology  optimization  problem  is  obtained  as  the  steady  state 
solution  of  the  coupled  system  (77),  i.e.  for  t  — >  oo  (or  T  sufficiently  large),  we  need  an 
adaptive  time  scheme  that  reduces  the  time  step  size  when  necessary  and  increases  it  as 
the  solution  approaches  the  steady  state.  We  employ  the  same  procedure  proposed  in 
[46]  for  the  Cahn-Hilliard  equation,  which  is  based  on  an  accuracy  criterion  and  reduces 
the  computational  costs  of  the  simulation  while  maintaing  an  adequate  level  of  accuracy. 

Very  recently,  a  new  second-order  accurate,  provably  unconditionally  stable,  time 
integration  algorithm  for  phase  Held  models  has  been  developed  in  [47].  This  would 
provide  a  viable  alternative  to  the  generalized-a  method. 


5.2.1  Time  step  scheme 

Let  us  subdivide  the  time  interval  [0,  T)  by  introducing  the  discrete  time  vector  {tn}"=0’ 
with  A tn  :=  tn. (-1  —  tn  the  width  of  the  time  interval  at  the  step  tn,  for  which  the  control 
variables  (82)  are  Pn+i  =  P(i„),  Pn+i  =  P (tn)  and  Un+i  =  U(t„).  Then,  if  we 
interpret  the  variables  Pn+i  and  P„+i  as  independent,  the  generalized-a  method  for 
the  TOft(D)  problem  at  time  step  tn+ \  reads  from  Eq.  (83): 


find  Pn_|_i,  Pn_|_i,  Pn+am  5  P n+a j  ?  Un+1 
Rp  P)i+aj-Tn+l  j  —  0; 

Ru  i  Pfi+a^  )  Un+1  j  0; 

Pn+1  =  Pri  +  At„P n  +  SAtn  ^P „_|_i  —  Pnj  , 

P?i+am  P n  +  ^Pn+1  P n  ^  5 

P n+a/  —  P n  4“  (Pn+1  P n)  5 


(84) 


with  Pn  and  Pra  given;  the  parameters  am,  ctf  and  S  £  K,  chosen  on  the  basis  of  stability 
and  accuracy  considerations,  define  a  specific  generalized-a  method.  As  described  in 
[54,  57]  we  select  am,  a/  and  S  as  follows: 


1  /  3  —  Poo\ 

2  U  +  Poo  )  ’ 


af  = 


1  +  Po 


6  —  ^  O:  f , 


(85) 


where  p <*,  G  [0, 1]  is  the  spectral  radius  of  the  amplification  matrix  at  At  -+  oo.  See 
[54,  57]  for  further  details. 

The  nonlinear  system  (84)  is  solved  for  each  time  step  t„+i,  for  n  =  0, . . . ,  nts  —  1,  by 
means  of  a  two  stage  predictor-multicorrector  algorithm,  for  which  the  control  variables 
at  the  time  step  tn+ 1  are  obtained  iteratively,  where  Pn+i,(j),  Pn+i,(»)  an<i  Un+1)(j), 
for  i  =  0, 1, . . . ,  imaxi  are  the  iterates  and  where  i  =  0  indicates  the  predictor.  At  the 
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predictor  stage  the  control  variables  are  initialized  as: 


71+1, (0) 


_ p 

5 


n+l,(0)  —  -trai 


Un+1,(0)  =  Un. 

At  the  multicorrector  stage  the  following  iteration  steps  are  repeated  for  i  =  1, 

1.  Update  the  control  variables  following  the  last  two  relations  of  Eq.  (84): 

Pn+am,(i)  =  P n  4“  |Pn+l,(i-l)  Pnj  5 

^n+o/,(i)  4“  af  (Pn+l,(i—  1)  Pn)  , 


Un+l,(t)  — 


2.  Assemble  the  residuals: 

Pp,(7)  *  Pp  (  Pji-{-am,(i)  7  Pn+a/ ,(i)  7  1  , 

) .  '  (88) 
Pu,(i)  *  Ru  (  Pn+am,(i)  1  Pfi-f  ap(i)  )  Un-j-qj)  J  • 

3.  If  the  following  stopping  criteria  on  the  relative  norms  of  the  residuals: 

l|Pp,(*)ll  ,7  ,  l|Pu,(i)||  , 

l|Pp,(0)ll  llKu,(0)ll 

are  satisfied  for  a  prescribed  tolerance  tolR,  set  the  control  variables  at  time  step 
ln+ 1  Pn+i  (*n+i,(i-i))  Pn+i  Pn+i,(i— l)  und  and  exit 

the  multicorrector  stage;  else  proceed  to  step  4. 

4.  Define  the  consistent  tangent  matrices  from  Eq.  (83): 


:=  Kee  ( 

/  . 

P 

A  n-\-cx.rr 

< 

Kpu>«) 

ii 

w 

c 

(  P 

l  A  n+ar 

K-u  p,(i) 

a 

II 

(  p 

l  A  n+ar 

:=KUU 

(  p 

\  A  n+Q;- 

n+ctm^i)  n-\-at f  ,(i)  n-\-l,(i)j  7 

n+am,(i)  t  ^n+a/  ,(i)  )  j  , 
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where,  by  using  Eq.  (84),  we  have: 


Knn  P 


(* 


pp  l  ^  n-\-0£m  5  ^  n-\-a.f  ?  Un_|_i 


K„n  P 


pu  1  ^  n+am  i  ^ 


Kuo  P 


(* 


up  l  ^ n-\-<y.m  7  ^ n+a/5  Un+1 


K, 


uu  l -^n+CKm  ? -^n+a/ 5  Un_|_i 


C^Pp  ^Pn+CKm  5  Pn+Qy  7  Un_)_l^  ^P 


n+a^ 


<9P 


<9P  n+Q;m 

<9R 

p  ^P n+a^  ?  P n+CK  j  i  U„+i)  ap 


n+l 


5RP  ( 


dP 


n+CK  / 


n+CK  f 

OPn+l 


Pn+ocm  5  Pn+aj  >  Un_|_i 


6  At, 
+af 


1  ^P  n-\-am 

<9RP  (P  n~\~otm  i  P n+CK  f  7  Un+lj 


dP 


n+c*  f 


dRn  P 


(* 


p  I  n+am  7  ^ n+a /  7  Un_|_i 


<9U 


n+l 


9RU( 


Pn+am5Pn+a/;Un+iy  9Pn+0 


<9P 


n-\-otn 


dP 


n+l 


+- 


SRU  (P 


n+am  >  fn+a /  j  Un+1  J  dPn+, 


dP 


(\rn  3Ru  ^Pn-(-om  ,  Pn-fa^  ,  Un+1 


n+a  / 


<9P 


n+l 


S  At, 
+af 


dP 


n+CK  77 


(  Pn+CKTTT,  7  Pn+a /  7  Un+1 


<9P 


n+CK  f 


dKu  (  Pn+CKm  7  Pn+a /  7  Un_|_i 


au 


n+l 


5.  Solve  the  following  linear  system  in  the  variables  AP„+1(j)  and  AU^+io): 

^^,(i)APn+l,(i)  “t“  Rp,(i) 

Rup,(i)^Pn+l,(z)  4”  R-uu,(?:)^Un+l,(i)  Ru,(i)* 


(91) 


(92) 


6.  Update  the  control  variables: 

Pn+l,(i)  =  Pn+l,(i-l)  +  ^Pn+l,(i), 

P  n+l,(i)  =Pn+l,(i-l)+AP„p:(i), 

Un+l,(t)  =  +  AUn+i^j). 

and  return  to  step  1. 
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5.2.2  Time  adaptivity 

We  consider  an  adaptive  scheme  similar  to  the  one  proposed  in  [46],  which  is  based 
on  the  comparison  of  the  solutions  obtained  with  the  generalized-a  method  and  the 
backward  Euler  method  [82].  The  backward  Euler  method  can  be  obtained  by  setting 
am  =  a/  =  S  =  lin  the  generalized-a  method. 

The  adaptive  time  scheme  starts,  for  each  time  step  tn+i,  n  =  0, . . .  ,nts  —  1,  with 
the  given  control  variables  Pn,  P„  and  Un,  and  a  given  time  step  A tn,  typically  that 
used  at  the  previous  time  step.  Then,  in  the  adaptive  algorithm  the  following  steps 
are  repeated  for  l  =  1  starting  with  A tn>(o)  =  A tn_i  (or,  if  n  =  0,  with 

Atnj(  o)  =  A  t0): 

1.  Compute  the  control  variables  Pn_|_ia_i),  Pn+i+j-i)  and  Ura+1(;_1)  with  the 
generalized-a  method  of  Sec.  5.2.1  for  A 

2.  Compute  the  control  variables  P^*^  P^+i  and  /^n  by  means  of 
the  backward  Euler  method  for  Atno_i). 

3.  If  the  generalized-a  or  the  backward  Euler  methods  are  not  converging  (i.e.  the 
predictor-multicorrector  algorithm  of  Sec.  5.2.1  is  not  convergent),  reduce  the  time 
step  by  means  of  a  safety  coefficient  xnc  G  (0, 1),  update  A tn^  as: 

Atn,(/)  =  Xnc  (94) 

and  return  to  step  1;  else  proceed  to  step  4. 


4.  Evaluate  the  relative  error  associated  to  A tnn_iy. 

IIPjj_i_r.a-.il  —  P^^i  n-.  IIU. 


eri+l,(£-l) 


X1-1)  ^n+l,(i-l) 

*71+1, (£  —  1)  II 


71+1, (£-1) 


-11  -u. 


BE 

77+1, (£-1)  I 


|U 


77+1, (£-1)1 


(95) 


5.  Update  the  time  step  size  according  to  the  following  formula: 


where: 


Xt7,(£— l)  :=  nun  <  xa 


tolA 


e7i+l,(£-l) 


(96) 

\ 1/2  1 

-J  ,  1  +  xgr  > 

(97) 

with  IoIa  a  prescribed  tolerance,  xa  €  (0, 1)  a  suitable  safety  coefficient  and 
Xgr  >  0  the  maximum  growth  rate  admitted. 


6.  If  e„+i  a_i)  >  IoIa ,  return  to  step  1.  Otherwise,  update  the  control  variables 
P77+1  i*77+i,(£— 1) 5  Pn+i  —  -Pt7+i,(£ — 1) 7  un+1  —  Un_j_r;p_ and  the  time  step 

Atn  =  A tj7,(i),  and  exit  the  loop. 

We  observe  that  on  the  basis  of  the  previous  algorithm,  one  of  following  three  situ¬ 
ations  occurs  after  step  6: 

•  if  e„+i,(£_i)  <  Xa^Ai  the  time  step  size  is  increased,  the  control  variables  of 
step  1  are  considered  valid  and  the  loop  ends; 

•  if  Xa^-oIa  A  e7i+i,(£-i)  <  tolA,  the  time  step  size  is  reduced,  the  control  variables 
of  step  1  are  still  considered  valid  and  loop  terminates; 
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•  if  e„+1(;_ i)  >  tolA ,  the  time  step  size  is  reduced,  the  control  variables  of  step  1 
are  invalid  and  the  steps  1-6  are  repeated. 

Since  the  steps  1  and  2  are  computationally  expensive,  the  occurrence  of  the  last  case 
should  be  minimum;  as  pointed  out  in  [46],  this  is  the  case  if  good  choices  for  the 
parameters  IoIa  and  xa  are  made. 

6  Selection  of  Parameters 

In  Sec.  4.3  we  provided  the  dimensionless  form  of  the  topology  optimization  prob¬ 
lem  (77),  which  we  proposed  to  solve  numerically  by  means  of  Isogeometric  Analysis 
and  the  generalized-a  method  in  Sec.  5.  However,  as  already  pointed  out,  the  solution 
of  the  TO(D)  problem  depends  on  the  dimensionless  parameters  D,  which  depend  on 
parameters  p,  £  V  introduced  in  Sec.  2.2.  In  this  section  we  address  the  issue  of  the 
choice  of  the  parameters  /igD  and  the  mesh  dependency  effect,  we  select  the  interpo¬ 
lation  function  g(p)  (13)  and  the  bulk  energy  V’s(p)  (16),  we  describe  the  continuation 
method  which  could  be  used  for  topology  optimization  in  the  phase  field  approach,  and 
we  discuss  the  choice  of  the  initial  solution  p0. 

6.1  The  choice  of  the  parameters  /teD  and  mesh  dependency 

The  dimensionless  parameters  D  (72)  depend  both  on  the  data  of  the  problem  and  the 
parameters  fi  =  (A,  7)  £  T>;  in  particular,  only  the  dimensionless  parameters  D2  and 
D3  and  the  representative  time  scale  To  are  related  to  p,  £  T>,  since  we  set  D 1  =  1. 

For  the  choice  of  the  parameter  A  we  adopt  similar  considerations  to  the  ones  made 
in  [46,  48].  The  Cahn-Hilliard  equation  converges  for  A  — >  0  (see  Eq.  (30))  to  a  thermo¬ 
dynamically  consistent  sharp-interface  model,  for  which  the  representative  length  scale 
of  the  interface  is  related  to  \/A.  Similarly,  for  the  generalized  Cahn-Hilliard  equation 
for  topology  optimization,  the  value  of  A  should  be  sufficiently  small  in  order  to  provide 
realistic  results  with  sharp  interfaces  between  the  material  and  the  void.  Moreover,  A  is 
responsible  for  the  number  of  holes  appearing  in  the  optimal  topology,  since  in  general 
it  is  associated  with  the  interface  energy  il>i(p)  of  Eq.  (17)  controlling  the  perimeter  of 
the  interfaces. 

From  a  numerical  point  of  view,  we  observe  that  the  computational  mesh  used  for 
the  spatial  approximation  should  be  fine  enough  to  capture  the  thin  layers  between  the 
material  and  void.  For  these  reasons,  following  from  the  dimensional  considerations  of 
Sec.  4.3,  we  assume  that  the  thickness  of  the  interfaces  depends  on  the  computational 
mesh  and  the  parameter  A  is  chosen  as  [46]: 


X  =  Xh2,  (98) 

where  h  is  the  characteristic  length  of  the  computational  mesh,  defined  as: 

h  :=  max  V^d,  (99) 

i=l,...,ne! 

with  Vi  the  area/volume  of  the  't-tli  element  of  the  mesh  composed  by  nei  elements. 
The  parameter  A  is  dimensionless  and  it  is  chosen  by  the  user.  Indeed,  it  is  difficult  to 
specify  at  this  point  what  value  A  should  assume  since  it  depends  in  general  on  the  data 
of  the  topology  optimization  problem;  for  example,  in  [46]  it  is  shown  that  for  the  Cahn- 
Hilliard  equation  A  should  depend  on  the  area/volume  V  covered  by  the  material  in  f l . 
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According  to  this  choice,  the  thickness  of  the  interfaces  is  roughly  equal  to  Vx  =  V Xh, 
with  shaper  interfaces  for  hner  meshes,  which  in  turn  allow  more  detailed  topologies  to 
be  obtained.  On  the  other  hand,  the  parameter  A  introduces  a  mesh  dependency  in  the 
topology  optimization  problem  through  D 2,  £>3  and  To;  this  represents  an  undesired 
issue  in  topology  optimization  since  the  optimal  distribution  of  the  material  in  the 
design  domain  O  changes  with  the  spatial  approximation  used.  However,  in  the  phase 
held  and  multiphase  approaches,  the  mesh  dependency  of  the  optimal  topology  could 
be  eliminated  by  using  a  fixed  value  of  h,  say  ho,  for  the  evaluation  of  the  parameter  A 
for  all  the  computational  meshes  having  h  <  ho- 

The  choice  of  the  parameter  7  is  important.  A  value  of  7  that  is  too  small  will 
cause  the  objective  functional  to  be  dominated  by  the  interface  and  bulk  energy  terms; 
if  7  is  too  large,  divergence  may  occur.  The  role  of  7  is  to  properly  balance  the  strain 
energy  function  ipE{p)  with  respect  to  the  bulk  and  interface  energy  functions,  V’b(p) 
and  tpi(p),  in  Eq.  (16).  However,  it  is  difficult  to  determine  a  priori  a  suitable  value  for  7 
by  means  of  dimensional  analysis;  this  is  due  to  the  fact  that  the  strain  energy  strongly 
depends  on  data  of  the  topology  optimization  problem,  which  are  not  represented  by  the 
dimensionless  parameters,  such  as  the  shape  of  the  design  domain  fl  and  the  direction 
of  the  body  and  surface  forces  f  and  h.  In  [114,  115],  the  values  differ  from  one  test  case 
to  another;  in  [23,  102]  the  choice  of  7  appears  to  be  arbitrary,  while  in  [20]  it  is  made 
by  trial  and  error.  In  order  to  ensure  a  proper  balance  of  the  strain  energy  with  respect 
to  the  other  energies  in  (22),  we  propose  to  decompose  the  dimensionless  parameter  7 
as: 

7  =  77 e,  (100) 

where  7  is  a  positive  parameter  chosen  by  the  user  and  7 E  is  computed  from  Eqs.  (72) 
and  (74)  as: 


7  e  — 


(A  A  (do)  +  A (do))  dfl 

[  D2'ipE(po,u(p0))dQ 


(101) 


with  po  the  initial  density,  ipE(po)  =  iPe(po-,  u(/°o))  from  Eq.  (69)  and,  in  order  to  make 
7 e  independent  of  A,  from  Eqs.  (72)  and  (98)  we  set 


A  :=  XD-2  = 


(102) 


Finally,  from  Eqs.  (98)  and  (101),  we  obtain  from  Eq.  (72)  that  the  dimensionless 
parameters  D2  and  Z)3  are: 

D 2  =  =^>.  As  =  77b  A  =  77  b==A  (103) 

A  hz  A/iA 

T4 

while  Tn  =  = — 2 —  since  we  have  set  D 1  =  1. 

Xh2M0 

At  this  point,  the  topology  optimization  problem  TO(D)  still  depends  on  the  arbi¬ 
trary  choices  made  for  the  parameters  A  and  7.  However,  numerical  tests  reveal  that 
these  dimensionless  parameters  vary  in  a  limited  range  of  values  for  different  topology 
optimization  problems;  indeed,  we  typically  select  A  €  [0.5,  6.0]  and  7  €  [0.5, 4.0].  On 
the  other  hand,  these  parameters  allow  the  user  to  modify  the  outcome  of  the  topology 
optimization  results.  In  particular,  the  number  of  holes  of  the  optimal  topology  can  be 
increased  or  decreased  by  decreasing  or  increasing  A,  respectively. 
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Figure  3:  Functions  f(p)  ( — )  and  fsiMp(p)  (-  -)  (left)  and  interpolation  functions  g{p) 
(— )  and  gsiAtp(p)  (-  -)  (right). 


6.2  The  choice  of  the  interpolation  function  g(p)  and  the  bulk 
energy  function  0 b(p ) 

As  anticipated  in  Sec.  2.1  the  interpolation  function  g(p)  is  typically  chosen  as  g(p)  =  pp 
(see  Eq.  (13)),  with  P  >  3  for  an  isotropic  material  with  vq  =  1/3.  We  consider  a  similar 
interpolation  rule  which  ensures  that  p  exceeds  a  minimum  value  and,  to  avoid 

numerical  issues,  the  condition  that  the  first  and  second  derivatives  of  such  interpolation 
function  are  zero  at  the  pure  phases  0  and  1.  In  particular,  we  consider: 


g{p)  =  f(p)p, 

with  f(p)  the  following  C2-continuous  function: 

(Pmin 

Pmin  (1  Pmin)^(p) 

1 


if  P<  Pi, 
if  pi  <  p  <  P6, 
if  P>  Pe, 


(104) 


(105) 


and  6(p)  a  C'2-continuous  third-degree  B-spline  in  p  obtained  with  the  open  knot  vector 
S  =  {{pi}t=1,P2,---,p5,{p6}t=i}  and  the  control  points  V  =  {0,  0,  0,  0, 1, 1, 1, 1};  see 
[80].  In  Fig.  3  we  compare  the  function  f(p)  with  the  typical  SIMP  function,  fsiMp(p), 
which  is  only  C-continuous  in  [0, 1],  where: 

{Pmin  if  P  <  Pmin i 

(106) 

P  if  P  A  Pmin  ? 


also,  we  compare  the  interpolation  function  g(p)  (104)  with  the  SIMP  function  say 
gsiMp(p)  =  fsiMp{p)P  for  P  =  3;  in  particular,  we  choose  pmin  =  0.1  and  S  = 
{{0.055}|=1, 0.065, 0.075, 0.925, 0.935,  {0.945}/=1}. 


29 


p 


Figure  4:  Bulk  energy  function  tpB^p). 


For  the  bulk  energy  function  iPb{p)  we  select  the  following  C°° -continuous  function 
in  p: 


V’s(p)  =  p2  (1  -  pf  +  A 


IQ —02p  _|_  ^q/32(p-1) 


(107) 


with  /3i,  /?2  €  Kq  .  This  choice  allows  naturally  steep  bounds  on  the  pure  phases  0  and 
1  but  avoids  any  singularities.  In  Fig.  4  we  plot  iPb(p)  for  p  in  [0, 1]  corresponding  to 
the  values  =  0.5  and  f32  =  50,  which  we  will  select  for  the  numerical  tests. 


6.3  The  continuation  method 

In  Remark  2.1  we  introduced  the  parameter  n  in  view  of  the  use  of  the  continuation 
method  to  solve  topology  optimization  problems  in  the  multiphase  approach  as  done 
for  example  in  [20,  23].  This  procedure  can  be  extended  to  the  case  of  the  topology 
optimization  problem  in  the  phase  field  approach  by  following  the  formulation  outlined 
in  Sec.  4  with  the  penalized  objective  functional  JK(p;/xfc)  (27).  The  dimensionless 
parameters  D\,  D2  and  D$  of  Eq.  (72)  are  modified  as: 

Dk  i  :=  —Di  DK  2  '■=  — 2^2,  Dk  3  :=  —D$,  (108) 

K  /-v  K 

for  which  DK  :=  (DKp,  -DK,2,  -D«,3,  £>4,  D5)  and  can  be  chosen  as: 

D*  :=  1kD2,  (109) 

for  some  7 K,  similar  to  Eq.  (103).  The  characteristic  time  is  TK  0  :=  —To,  corresponding 

K 

to  the  choice  DKp  =  1,  while  the  penalized  objective  functional  scales  with  the  quantity 

I  Ll~d , 

E  0  k\ 

We  propose  the  following  procedure  for  the  continuation  method: 

1.  set  the  data  for  the  topology  optimization  problem; 

2.  choose  an  initial  solution  po  such  that  /  po  dVL  =  V\ 

■In 

3.  select  the  parameters  A  and  7  of  Eqs.  (98)  and  (100),  respectively; 
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4.  choose  the  computational  mesh  for  the  spatial  approximation  and  compute  h  as 
in  Eq.  (103); 

5.  compute  the  parameter  A  from  Eq.  (98)  and  the  dimensionless  parameter  Di  from 
Eq.  (72); 

6.  select  a  discrete  set  K  of  values  for  the  parameter  k,  where  K  :=  with 

Ki  >  . . .  >  Km  >  0  and  M  the  number  of  continuation  levels; 

7.  repeat  the  following  operations  for  the  in  =  1 , ,M  continuation  levels: 

1.  if  m  =  1  set  p0,h,(i)  =  Po,h,  otherwise,  for  m  >  2,  set  p0,h,(m)  =  Ph,(m- 1)  for 
t  =  T,  where  Ph,(m- 1)  is  the  steady  state  solution  obtained  at  the  continuation 
level  to  —  1 ; 

2.  set  k  =  Km  and  compute  the  dimensionless  parameters  DK)(m)  from  Eqs.  (72), 

(108)  and  (109)  for  =  77 e,k,  where  7 e,k  is  computed  from  Eq.  (101)  by 

replacing  p0  with  p0,h,(m )  4;  also  compute  TKj0j(m); 

3.  solve  the  TO/,(DK](m))  problem  (81)  by  using  Isogeometric  Analysis  and  the 
generalized-ct  method  with  the  adaptive  scheme  outlined  in  Sec.  5  to  obtain 
the  steady  state  solution  Ph,(m)  for  t  =  T. 

The  choice  of  the  set  K  is  arbitrary  as  well  as  are  the  choices  of  the  parameters  A  and 
7.  As  a  general  indication,  if  the  number  of  continuation  levels  M  is  too  large,  high 
computational  cost  is  associated  with  the  solution  of  the  topology  optimization  problem 
with  this  approach.  On  the  other  hand,  if  M  is  too  small  or  the  ratios  Km-i / Km  1, 
then  convergence  may  not  occur. 

The  use  of  the  continuation  method  allows,  for  a  proper  selection  of  the  set  K,  to 
use  smaller  values  for  A  than  typically  used  when  the  TO^(D)  problem  (81)  is  solved  as 
a  single  simulation;  this  is  due  to  the  property  of  the  continuation  method  to  approach 
the  optimal  solution  as  a  sequence  of  intermediate  optimal  results. 

6.4  The  choice  of  the  initial  solution  p0 

In  the  phase  field  approach  the  topology  optimization  problem  allows  one  to  obtain  the 
optimal  solution  as  the  steady  state  of  the  phase  field  model  which  minimizes  a  penalized 
objective  functional  starting  from  an  initial  solution  po-  In  the  previous  sections  we  only 

required  that  po  satisfy  the  area/volume  constraint,  /  po  dCl  =  V.  It  follows  that,  even 

■In 

with  this  constraint  satisfied,  the  choice  of  p 0  is  arbitrary. 

Different  strategies  for  the  choice  of  po  have  been  considered  in  the  literature,  since 
the  initial  distribution  of  the  material  could  largely  influence  the  optimal  topology  in 
the  minimum  compliance  case  depending  on  the  approach  used.  In  general,  the  more 
Po  is  similar  to  the  optimal  solution,  the  more  rapid  will  be  the  convergence  of  the 
topology  optimization  method;  however,  this  situation  prefigures  an  a  priori  knowledge 
of  the  optimal  solution.  For  the  SIMP  approach  [12]  the  initial  solution  is  chosen  as 
Po  =  PVi  with  pv  '■=  V/|fi|,  even  if  po  often  is  allowed  to  violate  the  area/volume 
constraint,  which  is  later  restored  during  the  optimization  procedure.  In  the  level  set 
approach  [5,  35]  po  is  chosen  as  a  0-1  distribution  of  material  with  holes  such  that 
the  area/volume  constraints  is  satisfied;  the  initial  number  of  holes  affects  the  number 

4This  choice  allows  a  proper  balance  between  the  strain  energy  function  7 and  the  bulk  and 
interface  energies  ips  arid  7/  through  all  the  continuation  levels. 


31 


of  holes  in  the  optimal  topology  and  plays  a  crucial  role  in  the  definition  of  the  final 
topology.  However,  it  is  shown  in  [112]  by  means  of  numerical  examples  that  when  the 
capability  of  hole  nucleation  is  introduced  in  the  level  set  method,  the  dependence  of  the 
optimal  solution  on  the  initial  one  is  reduced.  Initial  solutions  in  the  0-1  configurations 
are  also  chosen  in  [102]  for  topology  optimization  problems  formulated  with  the  Cahn- 
Allen  equation,  which  shows  a  strong  dependence  of  the  optimal  solution  on  the  initial 
one.  For  the  phase  field  approach  of  [114,  115]  using  the  generalized  Cahn-Hilliard 
equation,  po  is  a  random  distribution  around  the  average  value  py. 

In  this  work,  we  typically  choose  po  =  Pv  ■  However,  the  continuation  method  of 
Sec.  6.2  can  be  effectively  used  to  provide  a  suitable  initial  solution  p0  without  any  a 
priori  knowledge  of  the  optimal  solution. 


7  Numerical  Tests 

In  this  section  we  solve  and  discuss  numerical  problems  in  two  and  three-dimensions;  we 

also  highlight  the  features  of  the  proposed  method  by  means  of  two-dimensional  tests. 

The  numerical  values  considered  for  the  numerical  simulations  and  the  implementation 

aspects  are  also  reported. 

7.1  Numerical  values  and  implementation  aspects 

For  all  the  simulations  reported  the  following  numerical  values  are  considered: 

•  for  the  tensor  Co  (see  Eq.  (1)),  we  choose  E0  =  200  •  109  J/m  and  vq  =  1/3  for 
the  plane  strain  two-dimensional  problems  and  the  three-dimensional  problems; 

•  for  the  representative  quantities  of  Sec.  4.3,  we  choose  Lq  =  1  m,  Tq  such  that 
Di  =  1,  Mq  =  1  m2/s,  ft®  =  200  •  106  Pa  (for  which  D4  =  1000)  and  /o  =  0  N/m3 
(no  body  force  is  considered  for  which  D$  =0); 

•  the  interpolation  g(p)  and  the  bulk  energy  V’b(p)  functions  are  chosen  as  in  Sec.  6.2 
with  P  =  3  and  the  mobility  is  assumed  constant,  M(p)  =  Mq. 

For  the  numerical  approximation  we  consider: 

•  for  the  Isogeometric  spatial  approximation  of  Sec.  5.1,  we  choose  a  NURBS  (or 
B-splines)  basis  of  degree  p  =  2  and  numerical  integration  is  performed  using 
3x3  Gauss  quadrature  (3  x  3  x  3  in  three-dimensional  problems);  this  represents 
a  very  conservative  approach.  The  study  of  more  efficient  quadrature  schemes  for 
NURBS  was  initiated  in  [56];  further  investigations  are  ongoing; 

•  for  the  time  step  scheme  of  Sec.  5.2.1  we  choose  p^  =  1/2  (for  which  am  = 
5/6,  cif  =  2/3  and  5  =  2/3  from  Eq.  (85)),  with  the  tolerance  for  the  stopping 
criterion  (89)  tola  =  10-4  and  the  maximum  number  of  steps  set  equal  to  imax  =  7; 

•  for  the  time  adaptivity  of  Sec.  5.2.2  we  select  the  initial  time  step  Afo  =  10~12 
and  the  final  time  T  =  102,  Xnc  ==  0-5  in  Eq.  (94)  to  reduce  the  time  step  in  case 
of  divergence  of  the  generalized-a  or  backward  Euler  methods,  and  \A  =  0.85, 
Xgr  =  1-2  and  tol^  =  10-3  in  Eq.  (97)  for  the  time  update  and  stopping  criterion; 
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•  the  numerical  solution  of  the  spatial  approximation  is  obtained  by  using  the  Bezier 
extraction  methods  presented  in  [17]  (Bezier  extraction  provides  a  localized  repre¬ 
sentation  of  the  globally  smooth  basis  that  can  be  implemented  in  shape  function 
routines  in  existing  finite  element  codes); 

•  the  problem  is  solved  using  a  parallel  C++  code  based  on  Trilinos  [53]  and  MPI; 

•  the  linear  system  (92)  is  solved  with  the  “monolithic”  approach  by  means  of 
the  GMRES  method  preconditioned  by  an  Algebraic  Multigrid  strategy  with 
Smoothed  Aggregation  [44];  the  dimension  of  the  Krilov  space  is  set  equal  to  700- 
1,000  and  the  stopping  criterion  is  based  on  the  relative  residual  with  tolerance 
equal  to  10~6  5. 

Remark  7.1  The  solution  of  the  linear  system  (92)  is  particularly  challenging  for  fine 
meshes,  especially  for  three-dimensional  problems;  this  is  due  to  the  nature  of  the  block 
matrices  involved  in  the  global  matrix.  Indeed,  the  matrix  K ppu)  radically  changes 
behavior  as  time  evolves:  when  A t  is  small,  far  from  the  steady  state,  the  mass  matrix 
dominates  over  the  fourth-order  stiffness  matrix  in  K pp,(i)  and  convergence  occurs  in 
a  relatively  small  number  of  GMRES  steps.  On  the  contrary,  when  At  is  large,  many 
more  steps  are  required.  The  block  matrix  Kuu  represents  the  stiffness  matrix  of  the 
linear  elastic  problem,  with  elastic  properties  depending  on  the  distribution  of  the  phase 
variable  p;  see  Eqs.  (1)  and  (10 f).  The  properties  of  such  a  matrix  abruptly  change  when 
passing  through  the  interfaces  and  GMRES  could  suffer  divergence  issues  when  sharp 
interfaces  are  generated,  especially  in  the  parallel  setting.  If  necessary,  a  way  to  partially 
overcome  this  inconvenience  consists  of  enforcing  the  development  of  the  interfaces  over 
a  sufficiently  high  number  of  mesh  elements  by  suitably  tuning  the  parameter  A  (98). 
Additionally,  the  off-diagonal  block  matrices  Kpu  m  and  KUP)^)  could  further  degrade 
the  conditioning  properties  by  rendering  the  full  matrix  “less  symmetric.  ”  We  believe 
that  an  ad  hoc  preconditioner  should  be  developed  to  take  into  account  all  these  features 
while  solving  the  linear  system  in  the  “monolithic”  approach;  we  notice  that  this  could  be 
required  even  if  the  system  (92)  is  eventually  solved  with  a  “staggered”  approach.  This 
was  not  pursued  in  the  present  work. 

On  the  basis  of  the  parameters  above  mentioned,  a  typical  two-dimensional  sim¬ 
ulation  with  a  mesh  with  320  x  160  elements  converges  in  about  nts  =  1,000-1,200 
time  steps;  time  adaptivity  converges  in  one  iteration  with  few  exceptions  for  a  limited 
number  of  time  steps  (typically  only  1  10).  In  each  time  step  the  generalized-a  method 
typically  converges  in  one  to  four  Newton  iterations  (likewise  for  the  backward  Euler 
method)  and  no,  or  a  few,  restarts  with  forced  reduction  of  the  time  step  are  required 
due  to  divergence.  For  the  same  simulation,  the  solution  of  the  linear  system  (92)  re¬ 
quires  a  number  of  GMRES  iterations  comprised  between  7  and  500  depending  on  the 
time  step  At  and  time  t.  In  general,  the  smaller  is  the  thickness  of  the  interface  for  a 
given  mesh  (related  to  the  parameter  A  of  Eq.  (98))  and/or  the  larger  the  parameter  7 
of  Eq.  (100),  the  slower  is  the  convergence  of  the  simulation  to  the  steady  state  solu¬ 
tion;  indeed,  in  these  cases,  we  could  incur  a  large  number  of  time  steps  and  time  step 

5  When  the  norm  of  the  residual  associated  to  the  solution  of  the  linear  system  is  below  the  10— 5 
threshold,  we  progressively  increase  the  tolerance  on  the  relative  residual.  This  situation  occurs  for 
“large”  At,  which  in  turn  occurs  in  proximity  of  the  steady  state  solution,  i.e.  when  the  solution  of  the 
problem  at  the  previous  time  step  yields  a  very  small  residual  in  the  linear  system  at  the  current  time 
step.  In  this  case,  a  fixed  tolerance  on  the  relative  residual  would  be  too  restrictive  and  would  lead  to 
an  unnecessary  large  number  of  GMRES  iterations. 
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Figure  5:  Test  1.  Design  domain  Q,  surface  force  h  and  displacement  constraints. 


adaptions,  and  a  large  number  of  Newton  iterations  in  the  generalized-a  and  backward 
Euler  methods  and,  consequently,  slow  convergence  of  the  GMRES  linear  solver. 

7.2  Two-dimensional  problems 

We  discuss  two-dimensional  topology  optimization  problems  in  order  to  highlight  the 
features  and  properties  of  the  method.  We  consider  design  domains  represented  by 
B-splines  and  NURBS  bases. 

7.2.1  Test  1 

The  topology  optimization  problem  for  Test  1  is  represented  in  Fig.  5.  We  consider  a 
rectangular  design  domain  fl  with  zero  displacement  boundary  conditions  on  the  left 
edge  and  the  vertical  traction  h  =  —  fi$n  acting  on  the  c^-part  of  the  right  edge.  We 
assume  L  =  2.00  m,  H  =  1.00  m,  d\  =  0.45  m  and  d2  =  0.10  m.  The  fraction  of  the 
volume  of  D  to  be  covered  by  the  material  is  T/|f2|  =  0.35.  The  design  domain  D  is 
represented  with  a  B-spline  basis  with  the  control  points  such  that  there  is  a  linear 
mapping  between  the  parametric  domain  and  Q. 

For  the  numerical  solution  of  the  problem  we  consider  a  mesh  with  320  x  160  elements, 
for  which  the  dimensionless  mesh  size  is  h  =  0.00625.  Also,  we  select  the  parameters 
A  =  6.0  and  7  =  1.0  for  the  definition  of  the  dimensionless  parameters  D2  and  D3  (see 
Eqs.  (98),  (100)  and  (103));  the  parameter  7 e  =  1.247  •  104  is  selected  from  Eq.  (101). 
The  dimensionless  parameters  are  D2  =  4.267  •  103  and  D3  =  5.322  •  107,  with  the 
characteristic  time  To  =  4.267  •  103.  We  indicate  this  problem  as  Test  1.1.  In  Fig.  6 
we  present  the  evolution  of  the  phase  (material  density)  variable  p  by  showing  it  at 
significant  time  steps;  the  initial  solution  is  set  as  po  =  Pv  =  V/|fi|.  The  color  scale  is 
from  blue  to  red  for  values  of  p  from  0  to  1.  As  we  can  deduce  from  Fig.  7 (left) ,  the 
distribution  of  the  material  is  driven  at  the  initial  steps  by  the  above-average  values 
of  the  strain  energy  function  1 pE]  in  this  case,  the  material  is  distributed  around  the 
r o  boundary  and  where  the  load  h  is  applied.  The  phase  separation  occurring  in 
these  areas  of  the  design  domain  is  due  to  the  generalized  Cahn-Hilliard  model  which 
is  mass/volume  conservative.  As  time  evolves,  the  separation  of  the  phases  is  mostly 
completed  and  the  structure  tends  to  simplify;  this  is  due  to  the  fact  that  the  bulk, 
interface  and  strain  energies  (i.e.,  Jb ,  Ji  and  Je,  respectively)  are  all  contributing  to 
the  minimization  of  the  objective  functional  J  (see  Eqs.  (18)  and  (24)).  The  steady  state 
represents  the  optimal  topology  in  terms  of  the  objective  functional  J,  which  reaches 
its  minimum  with  respect  to  the  initial  value  (for  p  =  po  at  t  =  0);  it  is  also  possible 
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t  =  1.644-  1CT6 


t  =  1.683-  1CT5 


t  =  6.528  •  10  4  (minimum  JE )  steady  state 

Figure  6:  Test  1.1.  Evolution  of  the  phase  (material  density)  variable  p  in  time  for  mesh 
size  320x160  with  A  =  6.0  and  7  =  1.0. 
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Figure  7:  Test  1.1.  Distribution  of  the  strain  energy  function  -!/Ae  (dimensionless)  in  the 
design  domain  D  (initial  and  at  the  steady  state);  logarithm  scale. 


to  see  from  Fig.  7(right)  that  the  average  value  of  the  strain  energy  function  ipE  is  also 
much  lower  at  the  steady  state.  However,  the  optimal  topology  at  the  steady  state 
does  not  necessarily  represents  the  stiffest  structure  with  the  specified  data;  indeed,  as 
is  occurring  in  this  case,  the  minimum  value  of  JE  is  obtained  at  an  intermediate  step 
and  before  the  steady  state  is  reached.  This  fact  is  highlighted  in  Fig.  8  where  the 
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Figure  8:  Test  1.1.  Normalized  energies  J  (black),  JE  (red),  JE  (blue)  and  Jj  (green) 
vs.  time  t  (left)  and  detail  (right);  the  minimum  value  of  JE  is  indicated  by  an  x. 


Figure  9:  Test  1.1.  Time  step  At  vs.  time  t  (dimensionless)  as  selected  by  the  adaptive 
time  scheme  (left)  and  detail  (right). 


behavior  of  the  objective  functional  J  and  the  energies  JE,  JE  and  J/  is  plotted  vs.  the 
dimensionless  time  t.  The  objective  functional  J  is  monotonically  decreasing  in  time  as 
expected  since  it  is  a  Liapunov  functional  (see  Eq.  (61)),  even  if  this  property  does  not 
hold  for  the  energies  JE,  JE  and  Jj.  The  large  drop  of  JE  is  due  to  the  phase  separation 
and  the  evolution  of  J/  to  the  creation  of  the  interfaces  and  the  following  simplification 
of  the  topology  during  the  evolution  in  time.  The  strain  energy  JE  is  also  subject  to 
a  large  and  desired  drop,  even  if  its  minimum  value  is  reached  at  an  intermediate  step. 
In  particular,  we  obtain  that  the  initial  value  of  the  objective  functional  J  decreases  by 
79.29%  at  the  steady  state,  while  JE  decreases  by  86.33%;  the  minimum  value  of  JE, 
representing  the  real  goal  of  the  topology  optimization  is  obtained  at  t  =  6. 528 TO”4  with 
a  drop  of  87.84%;  this  represents  a  significant  value  smaller  by  11.05%  with  respect  to 
the  one  reached  at  the  steady  state.  The  topologies  of  these  configurations  are  radically 
different;  consequently,  we  can  select  the  configuration  for  which  JE  is  minimum  as 
the  optimal  topology  from  a  design  point  of  view.  The  geometrical  information  can  be 
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A  =  6.0,  7  =  0.2 


A  =  12.0,  7  =  1.0 


minimum  Je 


minimum  Je 


steady  state 


steady  state 


Figure  10:  Test  1.2.  Phase  (material  density)  variable  p  for  mesh  size  320x160  with 
A  =  12.0  and  j  =  1.0  (left)  and  A  =  6.0  and  7  =  0.2  (right);  solutions  for  minimum 
value  of  Je  (top)  and  at  the  steady  state  (bottom).  The  contour  lines  represent  the 
topologies  obtained  in  Test  1.1  for  A  =  6.0  and  7  =  1.0. 


extracted  from  the  solution  by  the  contour  lines  corresponding  to  the  value  p  =  0.5. 
In  Fig.  9  we  show  the  evolution  of  the  time  step  A tn  vs.  time  t,  which  is  chosen 
adaptively  according  to  the  scheme  presented  in  Sec.  5.2.2.  We  observe  the  large  change 
in  the  value  of  A tn  from  10-12  to  18.8  occurring  in  a  relatively  small  number  of  time 
steps,  1, 198  in  this  case;  also,  we  observe  how  the  intermittent  nature  of  the  phase  field 
model  is  apparent  in  the  drop  of  the  value  of  A tn  required  at  some  times  steps.  As  a 
final  consideration  regarding  the  numerical  scheme,  we  observe  that  the  mass/volume 
constraint  is  adequately  satisfied  during  the  evolution  in  time,  since  the  relative  error 
with  respect  to  the  initial  value  of  the  volume  fraction  is  practically  negligible  and  below 
the  4.90  •  10-6%  threshold  for  all  the  time  steps. 

Finally,  we  discuss  the  results  of  the  topology  optimization  for  different  values  of  the 
parameters  A  and  7  with  respect  to  the  ones  previously  selected;  this  test  is  referred  to  as 
Test  1.2.  The  goal  is  to  show  that  a  large  value  of  A  leads  to  optimal  topologies  with  large 
interfaces  and  eventually  to  reduced  minimizations  of  the  strain  energy  Je-  Similarly, 
if  the  value  of  7  is  too  small,  the  phase  transition  leads  to  solutions  principally  driven 
by  the  Cahn-Hilliard  terms,  the  reduction  of  Je  would  be  limited  in  this  case  and  the 
topology  would  have  few  or  no  holes.  Ideally,  as  discussed  in  Sec.  6.1,  one  would  select 
the  smallest  possible  value  of  A  and  the  largest  of  7  which  would  lead  to  the  steady  state 
without  the  occurrence  of  divergence  issues  6.  In  Fig.  10(left)  we  present  the  topologies 
corresponding  to  the  minimum  value  of  Je  (top)  and  to  the  steady  state  (bottom)  for 
the  parameters  A  =  12.0  (twice  that  for  Test  1.1)  and  7  =  1-0;  the  contour  lines  for 

6  Divergence  issues  are  revealed  in  the  current  numerical  setting  by  increasing  values  in  time  of  the 
Liapunov  objective  functional  J  or  by  the  recursive  selection  of  excessively  small  values  of  the  adaptive 
time  step  A £n. 
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Figure  11:  Test  2.  Design  domain  fl,  surface  force  h  and  displacement  constraints. 


p  =  0.5  of  the  solutions  of  Test  1.1  are  highlighted.  We  can  observe  larger  interfaces 
with  respect  to  the  previous  case  (about  \/2  larger)  and,  even  if  the  configuration  at  the 
steady  state  is  very  similar  to  the  one  of  Test  1.1,  the  decrease  in  the  values  of  J  and  Je 
is  significantly  different.  Indeed,  the  drops  of  J  and  Je  at  the  steady  state  are  74.71% 
and  85.27%,  respectively  (for  Test  1.1  they  were  79.29%  and  86.33%);  at  its  minimum, 
the  drop  of  the  strain  energy  Je  is  86.61%,  a  value  larger  by  10.12%  with  respect  to  the 
one  obtained  in  Test  1.1.  Even  more  significant  differences  can  be  obtained  for  larger 
values  of  the  parameter  A.  In  Fig.  lO(right)  we  highlight  the  topologies  obtained  with 
the  parameters  A  =  6.0  and  7  =  0.2  (1/5  of  Test  1.1);  the  configurations  corresponding 
to  the  minimum  value  of  Je  and  at  the  steady  state  are  presented  together  with  the 
contour  lines  for  p  =  0.5  corresponding  to  Test  1.1.  The  thickness  of  the  interfaces 
between  the  phases  is  the  same  as  in  Test  1.1,  even  if  the  topologies  exhibit  significant 
changes.  Also,  the  relatively  “small”  value  of  7  leads  to  a  nonsymmetric  solution  at 
the  steady  state  since  the  symmetric  configuration  is  unstable  from  the  point  of  view  of 
the  phase  field  problem  (a  bifurcation  from  a  symmetric  configuration  occurs  during  the 
phase  transition).  The  decrease  of  J  and  Je  at  the  steady  state  are  79.62%  and  80.21%, 
respectively,  with  the  value  of  Je  being  34.35%  larger  than  in  Test  1.1;  similarly,  the 
drop  of  Je  at  its  minimum  is  83.05%,  39.39%  larger  than  in  Test  1.1.  Notice  that  in 
the  limit  case  7  =  0.0,  the  solution  coincides  with  the  Cahn-Hilliard  one,  which  in  this 
case  exhibits  a  topology  with  a  vertical  rectangle. 

7.2.2  Test  2 

For  Test  2,  we  consider  a  quarter  ring  with  inner  radius  Rin  =  1.00  m  and  outer  radius 
Rout  =  2.00  m  as  shown  in  Fig.  11.  The  surface  force  h  is  applied  at  the  top-left  tip  of 
the  design  domain  Q  along  a  boundary  segment  of  length  s  =  ir/lOm ;  no  displacement 
boundary  conditions  are  applied  at  the  bottom  of  Q.  The  volume  fraction  of  fi  to  be 
occupied  by  the  material  is  set  equal  to  V/|fi|  =  0.35.  We  represent  the  design  domain  fl 
by  means  of  a  NURBS  basis  of  degree  2  [80]  for  which,  by  using  Isogeometric  Analysis, 
the  exact  representation  of  the  geometry  is  maintained  through  the  analysis  and  the 
topology  optimization  procedure. 

We  solve  the  topology  optimization  problem  with  a  mesh  composed  by  320x160 
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t  =  9.933  •  10-6 


t  =  5.420  ■  10  5  (minimum  Je )  steady  state 

Figure  12:  Test  2.1.  Evolution  of  the  phase  (material  density)  variable  p  in  time  for 
mesh  size  320x160  with  A  =  5.0  and  7  =  1.5. 


elements  with  the  larger  number  along  the  circumferential  direction;  the  dimensionless 
mesh  dimension  is  h  =  9.818  •  10~3.  We  select  the  parameters  A  =  5.0  and  7  =  1.5,  for 
which  we  have  D2  =  2.075  ■  103  and  £>3  =  2.275  •  107  with  the  computed  7 e  =  7.309- 103. 
We  indicate  this  test  as  Test  2.1.  In  Fig.  12  we  present  the  evolution  in  time  of  the 
material  density  starting  from  the  initial  solution  po  =  pv  =  V/|fl|  toward  the  steady 
state;  intermediate  significant  solutions  are  presented,  including  the  one  corresponding 
to  the  minimum  value  of  the  strain  energy  Je-  We  observe  that,  even  if  such  solution 
does  not  exhibit  a  complete  phase  separation,  useful  topological  information  can  still  be 
obtained  and  the  configuration  can  be  used  for  further  design  investigation.  In  Fig.  13 
we  highlight  the  behavior  of  the  objective  functional  J  and  the  energies  Je ,  Jb  and  Jj 
vs.  the  time;  a  detailed  view  around  the  minimum  value  of  Je  is  also  presented.  The 
decrease  of  the  value  of  J  is  75.56%  at  the  steady  state,  while  for  Je  it  is  84.35%;  at 
its  minimum,  which  occurs  at  t  =  5.420  ■  10-5  (dimensionless),  the  decrease  of  Je  is 
86.25%,  12.12%  smaller  than  at  the  steady  state. 
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Figure  13:  Test  2.1.  Normalized  energies  J  (black),  JE  (red),  JE  (blue)  and  Jj  (green) 
vs.  time  t  (left)  and  detail  (right);  the  minimum  value  of  JE  is  indicated  by  an  x. 


We  now  discuss  the  mesh  dependency  effect  in  Test  2.2.  With  this  aim,  we  solve  the 
same  problem  discussed  previously  with  different  mesh  sizes  and  values  of  the  parameter 
A  affecting  the  dimensionless  parameters  D2  and  D3  (see  Eq.  (72)).  If  the  goal  of  the 
topology  optimization  is  to  obtain  sharp  interfaces  and  detailed  optimal  topologies,  fine 
meshes  need  to  be  used  since,  from  Eq.  (98),  the  parameter  A  is  selected  as  A  =  Ah2,  with 
h  indicative  of  the  mesh  size.  However,  as  already  mentioned  in  Sec.  6.1,  this  introduces 
a  mesh  dependency  effect  on  the  solution,  since  the  dimensionless  parameters  D2  and 
D3  vary  with  the  mesh.  On  the  other  hand,  the  mesh  dependency  issue  can  be  quickly 
eliminated  in  the  phase  field  approach  by  using  a  fixed  value  of  h  =  ho  for  all  the 
meshes  having  h  <  ho-  In  Fig.  14(right)  we  show  the  optimal  topologies  obtained  at 
the  steady  states  for  different  mesh  sizes  (80x40,  160x80  and  320x160)  for  A  =  5.0  and 
7  =  1.5;  once  ho  is  set  to  be  the  representative  dimensionless  size  of  the  mesh  80x40 
(h0  =  0.03927),  we  select  h  =  ho,  ho/2,  ho/4  for  the  three  meshes,  respectively.  In  this 
case,  we  observe  that,  not  only  does  the  thickness  of  the  interfaces  change  from  one  mesh 
to  the  other,  but  also  the  optimal  topologies  significantly  vary.  Conversely,  if  we  assume 
a  constant  value  of  h  =  ho  for  all  the  meshes,  the  thickness  of  the  interfaces  between  the 
phases  and  the  optimal  topologies  remain  the  same;  see  Fig.  14(left)  where  we  considered 
A  =  1.0  and  7  =  2.0  to  show  shaper  interfaces  (the  same  result  of  Fig.  14(top-right) 
would  have  been  obtained  for  all  three  meshes  with  A  =  5.0  and  7  =  1-5).  These  facts 
are  better  highlighted  in  Fig.  15  in  terms  of  the  behavior  of  the  objective  functional  J 
and  strain  energy  JE  for  the  three  meshes  with  and  without  a  fixed  value  of  h.  As  we 
can  observe  in  Fig.  15(left)  the  choice  of  a  fixed  h  =  ho  leads  to  a  good  match  between 
the  energies  for  all  the  meshes,  with  only  minor  differences;  at  the  steady  state  the 
maximum  discrepancy  on  J  with  respect  to  the  finer  mesh  is  0.03630%,  and  0.4021% 
for  JE.  On  the  contrary,  the  mesh  dependency  effect  is  clearly  visible  in  Fig.  15(right) 
with  the  mesh  dependent  h  =  ho,  ho/2,  ho/ A.  Additionally,  we  observe  that  the  finer 
meshes  allow  a  more  significant  reduction  of  the  objective  functional  J  and  strain  energy 
J E  with  respect  to  the  coarse  one,  due  to  the  ability  to  deal  with  sharper  interfaces. 
In  particular,  the  reduction  of  J  at  the  steady  state  is  51.56%,  67.93%  and  75.56%  for 
the  three  meshes,  while  for  JE  it  is  77.73%,  83.63%  and  84.35%;  at  their  minimum  the 
reductions  of  JE  are  78.46%,  83.79%  and  86.25%,  with  the  maximum  stiffness  obtained 
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mesh  80x40,  h  =  ho 


mesh  80x40,  h  =  ho 


mesh  320x160,  h  =  ho  mesh  320x160,  h  =  ho /A 

Figure  14:  Test  2.2.  Steady  states  of  the  phase  (material  density)  variable  p  for  the 
mesh  sizes  80x40,  160x80  and  320x160  with  fixed  h  =  ho  (A  =  1.0,  7  =  2.0)  (left)  and 
with  mesh  dependent  h  =  ho,  h0/ 2  and  h0/A  (A  =  5.0,  7  =  1.5)  (right);  h0  =  0.03925. 
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Figure  15:  Test  2.2.  Normalized  energies  J  (black)  and  Je  (red)  vs.  time  t  for  the 

mesh  sizes  80x40  (-  •),  160x80  ( - )  and  320x160  ( — )  for  fixed  h  =  h.0  (A  =  1.0, 

7  =  2.0)  (left)  and  mesh  dependent  h  =  ho,  ho/2  and  h0/4  (A  =  5.0,  7  =  1.5)  (right). 


with  the  finer  mesh. 

Finally,  we  solve  the  topology  optimization  problem  by  means  of  the  continuation 
method  described  in  Sec.  6.3;  the  goal  is  to  show  that  this  approach  can  be  conveniently 
used  to  generate  solutions  with  sharp  interfaces  even  if  a  coarse  mesh  is  used.  We 
refer  to  this  problem  as  Test  2.3.  With  this  aim,  we  consider  a  mesh  of  size  80x40 
(dimensionless  h  =  0.03927  from  Eq.  (99))  with  a  two-level  continuation  procedure  for 
which  the  parameters  n  are  chosen  as  k  €  K  =  {4.0, 1.0};  additionally,  we  select  A  =  1.0 
and  7  =  4.0  (see  Eqs.  (98)  and  (100)).  The  resulting  dimensionless  parameters  (108)  are: 
.Dk, 2  =  4.053- 101  and  DKt3  =  4.742- 106,  with  7 e,k  =  7.313 - 103,  for  continuation  level  1, 
and  Dk^ 2  =  6.485  •  102  and  3  =  7.544  •  107,  with  7 e,k  =  2.907  ■  104,  for  level  2.  In 
Fig.  16  we  present  the  result  of  the  topology  optimization  with  the  two-level  continuation 
method  in  which  we  highlight  the  steady  states  as  well  as  the  solution  at  significant 
time  steps  for  levels  1  and  2;  the  optimal  topology  is  represented  by  the  steady  state  of 
level  2.  The  evolution  of  the  normalized  energy  J  in  time  t  (dimensionless)  is  presented 
in  Fig.  17  for  the  two  levels.  The  total  and  elastic  energies  decrease  72.16%  and  82.09% 
for  continuation  level  1,  and  27.68%  and  35.64%  for  level  2;  for  the  whole  procedure, 
the  decrease  of  J  and  Je  with  respect  to  the  values  corresponding  to  the  initial  solution 
Po  =  Pv  is  79.87%  and  88.47%,  respectively.  For  comparison,  we  observe  that  divergence 
issues  appear  for  the  topology  optimization  of  the  current  test  without  the  continuation 
method  for  A  =  1.0  and  7  =  4.0  (the  same  values  used  for  the  continuation  method  with 
D2  =  DKi 2  for  k  =  1.0).  It  follows  that  in  order  to  maintain  the  same  thickness  of  the 
interfaces  of  the  solution  as  for  the  continuation  method,  the  value  of  7  would  need  to 
be  lowered  (the  value  7  =  2.0  would  suffice  for  the  current  mesh),  even  if  the  reduction 
of  Je  would  be  smaller  than  with  the  continuation  method. 


42 


Figure  16:  Test  2.3.  Evolution  of  the  phase  (material  density)  variable  p  in  time  t  for 
mesh  size  80x40  with  the  continuation  method  for  K  =  {4.0, 1.0},  A  =  1.0  and  7  =  4.0; 
continuation  levels  1  (top)  and  2  (bottom)  with  steady  states  (right). 
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Level  1  Level  2 

Figure  17:  Test  2.3.  Normalized  energies  J  (black),  Je  (red),  Jg  (blue)  and  Jj  (green) 
vs.  time  t  for  the  continuation  method  at  levels  1  (left)  and  2  (right). 
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Figure  18:  Test  3.  Design  domain  fl,  surface  force  h  and  displacement  constraints. 


7.2.3  Test  3 

For  this  test  case,  we  consider  the  topology  optimization  problem  depicted  in  Fig.  18. 
The  surface  forces,  h,  are  applied  over  segments  of  the  boundary  of  length  d3.  We 
enforce  zero  displacement  condition  on  a  segment  of  the  boundary  on  the  left  end  of  the 
bottom  edge  and  we  fix  the  vertical  displacement  over  a  segment  on  the  right  end  of 
the  bottom  edge  7.  We  assume  L  =  4.00  m,  H  =  1.00  m,  d\  =  0.125  m,  c?2  =  0.9375  m, 
d3  =  0.125  m  and  d 4  =  0.875  m;  the  volume  fraction  to  be  occupied  by  the  material 
is  V/|f2|  =  0.40.  The  design  domain  D  is  represented  by  means  of  a  B-splinc  basis  of 
degree  2  with  a  linear  mapping. 

With  this  test  problem,  we  discuss  the  solution  of  topology  optimization  problems 
in  the  presence  of  local  peak  values  of  the  strain  energy  function  i/je  inside  the  design 

'  In  the  framework  of  Isogeometric  analysis,  we  impose  the  displacement  constraints  on  the  degrees 
of  freedom  corresponding  to  the  control  points  on  the  Dirichlet  boundary  segments  Tg. 
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Figure  19:  Test  3.  Distribution  of  the  strain  energy  function  ipE  (dimensionless)  in  the 
design  domain  D  at  t  =  0  for  po  =  py  with  mesh  size  128x80;  logarithm  scale. 


Level  1,  steady  state 


Level  2,  steady  state 

Figure  20:  Test  3.1.  Evolution  of  the  phase  (material  density)  variable  p  in  time  t  with 
the  continuation  method  for  K  =  {8.0, 1.0},  A  =  2.0  and  7  =  1.0;  continuation  levels  1 
(top)  and  2  (bottom)  at  the  steady  states. 


domain.  An  example  is  highlighted  in  Fig.  19,  where  the  distribution  of  1 pE  is  shown 
for  this  test  problem  at  the  initial  step  (po  =  pv  =  V/|fi|).  In  such  cases,  the  phase 
separation  is  rapid  and  locally  driven  by  such  peak  values  and  convergence  may  not  occur 
if  the  values  of  the  parameter  7  are  too  “large”  and/or  those  of  the  parameter  A  are  too 
“small.”  However,  since  the  goal  is  to  obtain  a  significant  reduction  of  the  strain  energy 
Je  with  sharp  interfaces,  the  choice  of  the  parameters  cannot  be  excessively  limited  by 
such  issues.  In  practice,  “small”  values  of  7  and  “large”  values  of  A  should  be  chosen 
at  the  early  stages  of  the  phase  transition  and  modified  during  its  evolution.  However, 
we  realize  that  this  procedure  needs  to  be  calibrated  for  each  topology  optimization 
problem.  In  order  to  overcome  this  deficiency,  we  propose  two  approaches.  The  first 
one  consists  in  using  the  continuation  method,  while  the  second  one  considers  an  initial 
phase  (material)  distribution  po  in  which  “bubbles”  of  material  are  located  ab  initio 
in  correspondence  with  peak  values  of  the  strain  energy  function  %pE-  The  mesh  size 
chosen  for  the  simulations  is  512x128  and  h  =  0.007813. 

In  Fig.  20  we  present  the  solution  of  the  topology  optimization  problem  using  the 
continuation  method;  the  two  steady  state  solutions,  corresponding  to  the  values  k  = 
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Level  2 


Figure  21:  Test  3.1.  Normalized  energies  J  (black),  Je  (red),  Jg  (blue)  and  Jj  (green) 
vs.  time  t  for  the  continuation  method  at  levels  1  (left)  and  2  (right). 


{8.0, 1.0},  are  presented.  For  this  simulation,  which  we  indicate  as  Test  3.1,  we  choose 
A  =  2.0  and  7  =  1.0.  We  obtain  the  following  dimensionless  parameters  from  Eq.  (108): 
D K  2  =  1-280  •  102  and  DK} 3  =  1.074  •  107,  with  7 e,k  =  1-049  •  104,  for  continuation 
level  1,  and  Dk2  =  8.192  •  103  and  3  =  1.989  •  10s ,  with  7 e,k  =  2.428  •  104,  for  level  2. 
In  Fig.  21  we  present  the  normalized  energies  vs.  the  dimensionless  time  t.  The  overall 
decrease  of  the  total  and  strain  energies  with  respect  to  the  initial  solution  p 0  =  pv  is 
81.31%  and  85.14%,  respectively  (the  decrease  of  J  and  Je  at  the  continuation  level  1 
is  72.49%  and  81.54%,  while  at  level  2  it  is  38.07%  and  19.48%).  The  minimum  value 
of  Je  is  obtained  at  the  steady  state  of  continuation  level  2. 

In  Fig.  22  (bottom)  we  present  the  steady  state  of  the  phase  (material)  variable  ob¬ 
tained  by  solving  the  topology  optimization  problem  starting  from  the  initial  solution  p3 
depicted  in  Fig.  22(top);  notice  the  “bubbles”  of  material  distributed  in  correspondence 
of  the  peak  values  of  the  strain  energy  function  i/je  (see  Fig.  19).  This  test  is  referred  as 
Test  3.2.  The  volume  fraction  of  this  case  is  V/|f2|  =  0.4067  due  to  the  presence  of  the 
initial  “bubbles”;  also,  we  choose  A  =  6.00  and  7  =  1.00.  The  dimensionless  parameters 
of  Eq.  (72)  are  D2  =  2.731  •  103  and  D3  =  3.452  •  107  with  7^  =  1.264  •  104.  In  Fig.  23 
we  present  the  behavior  of  the  normalized  energies  with  respect  the  dimensionless  time 
t.  The  overall  reductions  of  J  and  Je  at  the  steady  state  are  77.58%  and  80.85%,  re¬ 
spectively.  The  maximum  reduction  of  Je  (82.03%)  is  obtained  at  t  =  1.818  •  10-4;  the 
corresponding  solution  is  shown  in  Fig.  22  together  with  another  significant  topology 
obtained  during  the  phase  transition  (at  t  =  9.217  •  10~5  the  reduction  of  Je  is  81.74%, 
4.650%  smaller  than  at  the  steady  state),  which  can  be  selected  for  further  design  and 
analysis  investigation. 
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t  =  0 


t  =  9.217-  10”5 


t  =  1.818  •  10  4  (minimum  Je) 


steady  state 

Figure  22:  Test  3.2.  Evolution  of  the  phase  (material  density)  variable  p  in  time  t  with 
A  =  6.0  and  7  =  1.0;  prescribed  initial  solution  p0  with  “bubbles”  (top)  and  steady 
state  (bottom). 
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Figure  23:  Test  3.2.  Normalized  energies  J  (black),  JE  (red),  JE  (blue)  and  Jj  (green) 
vs.  time  t  (left)  and  detail  (right);  the  minimum  value  of  JE  is  indicated  by  an  x. 


Figure  24:  Test  4.  Design  domain  f 2,  surface  force  h  and  displacement  constraints. 


7.3  Three-dimensional  problems 

In  this  section,  we  solve  topology  optimization  problems  defined  in  three-dimensional 
design  domains  fi. 

7.3.1  Test  4 

We  consider  the  topology  optimization  problem  represented  in  Fig.  24  in  which  the 
design  domain  D  is  a  solid  beam  of  size  HxHxL ,  the  surface  force  h  =  —  Hoy  is  applied 
on  a  subdomain  of  the  front  face  and  zero  displacements  are  imposed  on  the  back  face 
(the  plane  z  =  0).  We  assume  L  =  4.00m,  H  =  1.00m,  d  =  0.375 m  and  s  =  0.250m; 
the  volume  fraction  to  be  occupied  by  the  material  is  set  equal  to  V/|ST|  =  0.35.  By 
using  symmetry  properties,  only  half  of  the  domain  f 1  (whose  size  is  H/2xHxL)  is 
considered  for  the  computations.  The  design  domain  is  represented  by  a  B-splines  basis 
of  degree  2  composed  of  10x20x80  elements.  To  enforce  the  symmetry  of  the  solution 


48 


t  =  4.401  •  10"5 


t  =  2.577  •  10~4 


t  =  5.760  •  10  3  (minimum  Je) 


t  =  1.940  •  10”1 


t  =  1.946  •  10  1  steady  state 

Figure  25:  Test  4.1.  Evolution  of  the  phase  (material  density)  variable  p  in  time  t  for 
mesh  size  10  x  20  x  80  with  A  =  1.0  and  7  =  1.5;  the  volume  to  be  occupied  by  the 
material  (for  p  >  0.5)  is  displayed;  the  color  gray  indicates  p  =  0.5. 


with  respect  to  the  plane  y  =  H/2  and  prevent  the  detection  of  a  nonsymmetric  local 
minimum,  we  impose  the  same  values  of  the  phase  variable  on  the  planes  y  =  0  and 
y  =  H. 

We  solve  the  problem  by  considering  A  =  1.0  and  7  =  1.5  with  the  dimensionless 
size  of  the  mesh  h  =  0.05  (Test  4.1);  the  computed  parameter  7 e  takes  the  value  7 e  = 
2.040-104.  The  corresponding  values  of  the  dimensionless  parameters  are  D2  =  4.000-102 
and  D3  =  1.224  •  107  with  the  characteristic  time  To  =  4.000  •  102.  In  Fig.  25  we  depict 
the  phase  variable  p  at  significant  time  steps  (dimensionless),  including  the  steady  state 
solution  and  the  solution  corresponding  to  the  minimum  value  of  Je;  the  regions  of  the 
design  domain  Q  occupied  by  material,  which  are  obtained  for  p  >  0.5,  are  highlighted. 
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minimum  Je 


steady  state 


Figure  26:  Test  4.1.  Optimal  topologies  corresponding  to  the  minimum  value  of  Je  (left) 
and  to  the  steady  state  (right)  obtained  for  A  =  1.0  and  7  =  1.5. 


Figure  27:  Test  4.1.  Normalized  energies  J  (black),  JE  (red),  Jb  (blue)  and  Jj  (green) 
vs.  time  t  (left)  and  detail  (right);  the  minimum  value  of  Je  is  indicated  by  an  x. 


The  capability  of  the  method  to  handle  hole  nucleation  is  highlighted  in  Fig.  25  for 
the  solutions  at  the  time  steps  t  =  1.940  ■  10-1  and  t  =  1.946  •  10”1,  in  which  a  hole 
(specific  volume  of  the  design  domain  with  p  <  0.5)  is  generated.  In  Fig.  26(left)  we 
show  the  optimal  configuration  for  which  the  value  of  Je  is  minimum;  in  Fig.  26(right) 
the  configuration  at  the  steady  state  is  shown  for  comparison.  The  evolution  of  the 
objective  functional  J  and  the  energies  Je,  Jb  and  J/  vs.  the  dimensionless  time  are 
shown  in  Fig.  27(left);  in  Fig.  27(right)  a  detail  around  the  minimum  value  of  Je  is 
presented.  At  the  steady  state,  the  drops  of  the  values  of  J  and  Je  with  respect  to  the 
initial  solution  po  =  Pv  are  75.78%  and  87.38%,  respectively.  The  minimum  compliance 
is  obtained  at  t  =  5.760  •  10~3  with  an  88.66%  decrease;  at  its  minimum,  the  value  of 
Je  is  10.13%  smaller  than  at  the  steady  state.  In  general,  all  the  topologies  obtained  in 
the  range  t  =  1.200- 10-3  -  10-2  could  be  eventually  considered  for  further  analysis  and 
design  investigation  since  the  values  of  Je  are  close  to  the  minimum  (less  than  2.00% 
difference) . 

For  comparison,  we  solve  the  same  problem  with  different  parameters  A  and  7. 
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Figure  28:  Test  4.2.  Optimal  topology  corresponding  to  the  minimum  value  of  Je  for 
A  =  0.75  and  7  =  2.0. 


For  this  test,  which  we  indicate  as  Test  4.2,  we  consider  A  =  0.75  and  7  =  2.0;  the 
dimensionless  parameters  are  D2  =  5.333  •  102  and  O3  =  2.176  •  10' .  Due  to  the  larger 
value  of  7  and  the  smaller  value  of  A,  we  expect  a  more  significant  reduction  of  the 
strain  energy  Je  with  slightly  thinner  interfaces.  In  Fig.  28  we  present  the  optimal 
configuration  corresponding  to  the  minimum  value  of  Je,  which  is  14.16%  less  than  the 
one  obtained  for  Test  4.1. 
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Figure  29:  Test  5.  Design  domain  12,  surface  force  h  and  displacement  constraints. 


7.3.2  Test  5 

This  problem  is  defined  in  a  design  domain  12  represented  by  an  hemispherical  thin  shell 
as  shown  in  Fig.  29.  The  surface  force  h  =  and  the  displacement  constraints 

are  also  shown.  We  assume  R  =  0.250 m,  H  =  1.00m,  d  =  0.0314 m  (arc  length), 
s  =  0.200  m  (arc  length)  and  thickness  a  =  0.0200  m;  the  inner  radius  of  the  shell  is 
constant  and  equal  to  H  for  all  the  planes  through  the  axis  y.  The  thin  shell  is  described 
by  a  single  NURBS  patch  with  basis  of  degree  p  =  2  [80];  by  virtue  of  the  symmetry 
properties,  only  a  quarter  of  the  shell  is  considered  as  the  domain  for  the  computations. 
The  volume  fraction  to  be  occupied  by  the  material  is  V/|12|  =  0.35. 

The  shell  is  modelled  as  a  three-dimensional  linear  elasticity  problem.  For  the 
computations  we  consider  a  mesh  composed  of  1282xl  elements.  In  this  manner,  since 
we  consider  a  basis  of  degree  p  =  2,  the  bending  modes  are  properly  taken  into  account 
by  the  three-dimensional  linear  elastic  model  [13];  indeed,  three  control  points  exist 
through  the  thickness.  An  additional  consequence  of  this  choice  is  that  the  distribution 
of  the  phase  variable  p  is  constant  throughout  the  thickness  due  to  the  conditions  Vp-n  = 
0  on  the  opposite  faces.  In  order  to  enforce  the  symmetry  of  the  solution  in  the  quarter 
of  the  shell,  we  impose  the  same  values  of  the  phase  variable  at  the  x  =  0  and  z  =  0 
planes. 
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Level  2,  t  =  1.596  •  10  6  Level  2,  steady  state 


Figure  30:  Test  5.  Evolution  of  the  phase  (material  density)  variable  p  in  time  t  for 
mesh  size  1282  x  1  with  the  continuation  method  for  IK  =  {5.0, 1.0},  A  =  1.0  and  7  =  2.5; 
continuation  levels  1  (top)  and  2  (bottom)  with  steady  states  (right). 


We  solve  the  topology  optimization  problem  by  means  of  the  continuation  method 
with  two  levels;  see  Sec.  6.3.  In  particular,  we  assume  k  €  1  =  {5.0, 1.0},  A  =  1.0 
and  7  =  2.5;  the  dimensionless  element  size  is  chosen  as  h  =  0.01535.  The  resulting 
dimensionless  parameters  are:  Dk  2  =  1.698  ■  102,  UKj3  =  9.770  •  105,  7 e,k  =  4.604  •  102 
for  continuation  level  1,  and  DK  2  =  4.244  •  103,  Dk  3  =  1.404  •  10' ,  7 e,k  =  1-323  ■  103  for 
level  2.  In  Fig.  30  we  present  the  evolution  of  the  material  density  in  time  throughout 
the  two  continuation  levels.  The  final  and  optimal  topology  is  the  steady  state  of  the 
second  continuation  level.  In  Fig.  31  we  show  the  optimal  configuration  on  the  whole 
design  domain  obtained  by  placing  the  material  where  p  >  0.5.  In  Fig.  32  we  show 
the  evolution  of  the  objective  functional  J  and  the  energies  Je,  Jb  and  Jj  vs.  the 
dimensionless  time  for  both  continuation  levels.  Specifically,  we  obtain  drops  of  71.83% 
and  79.50%  for  J  and  Je  for  level  1,  and  27.17%  and  28.83%  for  level  2,  respectively. 
The  overall  decreases  of  J  and  Je  throughout  the  whole  continuation  level  procedure 
are  79.48%  and  85.41%,  respectively;  the  minimum  compliance  is  obtained  at  the  steady 
state  of  continuation  level  2. 
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Figure  31:  Test  5.  Optimal  topology  obtained  at  the  steady  state  of  continuation  level  2. 


Figure  32:  Test  5.  Normalized  energies  J  (black),  Je  (red),  Jb  (blue)  and  Jj  (green) 
vs.  time  t  for  the  continuation  method  at  levels  1  (left)  and  2  (right). 


8  Conclusions 

In  this  work  we  solved  minimum  compliance  topology  optimization  problems  with  a 
phase  field  model  based  on  the  generalized  Cahn-Hilliard  equation.  With  this  ap¬ 
proach,  the  interfaces  between  the  phases,  material  and  void,  are  represented  by  sharp, 
but  smooth,  layers  and  the  optimal  solution  is  obtained  as  the  steady  state  of  the  phase 
transition  problem,  eliminating  the  need  of  an  optimizer.  The  ability  to  deal  with  topo¬ 
logical  changes  and  hole  nucleation  is  naturally  embedded  in  the  model  as  well  as  the 
constraint  on  the  total  amount  of  material  to  be  distributed  in  the  design  domain.  With 
this  formulation,  the  topology  optimization  problem  is  completely  defined  at  the  con¬ 
tinuous  level  and  the  optimal  solution  depends  on  dimensionless  parameters  specified 
by  the  user.  Additionally,  the  mesh  dependency  effect,  which  typically  affects  the  opti¬ 
mal  topologies,  can  be  eliminated  by  a  suitable  choice  of  the  parameters  controlling  the 
thickness  of  the  interfaces  and  the  number  of  holes  in  the  topology.  The  continuation 
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method,  a  multilevel  optimization  strategy  often  used  for  the  solution  of  the  topology 
optimization  problems,  is  extended  to  the  phase  field  approach. 

For  the  numerical  approximation  we  used  Isogeometric  Analysis,  which  is  particularly 
suitable  for  phase  field  problems  and  allows  exact  CAD  geometry  to  be  used  to  describe 
the  design  domain  and  to  also  be  used  in  the  optimization  procedure.  For  the  time 
approximation  we  used  the  generalized-a  method  in  combination  with  a  time-adaptive 
scheme  which  allowed  to  efficiently  capture  the  fast  and  intermitted  variations  in  time 
typically  occuring  in  the  phase  field  model. 

We  solved  both  two  and  three-dimensional  problems  to  illustrate  the  validity  of  the 
approach,  which  we  believe  is  a  promising  one  for  topology  optimization  problems. 
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